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Title:  A  Lagrangian  Fractional  Step  Method  for  the  Incompjressible  Navier-Stokes 
Eqxiatiom 

Aathor:  Christoph  BOrgers 

AdvlBor:  Charles  S.  Pesldn 

Abstract:  We  develop  a  modification  of  Peskin's  Lagrangian  fractional  step  method  for 
the  incompresaibk  Navier-Stokes  equations.  This  now  method  is  substantially  more  effi- 
cient than  the  one  originally  proposed  \  .r«kin.  On  a  grid  with  N  points,  the  work  per 
time  step  is  proportional  to  NlogN,  This  gain  in  efficiency  is  accomplished  by  modifying 
the  splitting  and  by  using  a  multigrid  method  for  the  solution  of  the  resulting  systems  of 
equations. 

The  method  uses  finiie  difference  operators  constructed  with  the  aid  of  Voronoi 
diagrams.  We  have  imi  <!etri.ented  it  on  a  periodic  domain  in  the  plane.  We  describe  an 
efficient  algorithm  for  the  numf;rical  con:  ituction  of  periodic  Voronoi  diagrams  in  the 
plane.  We  believe  that  this  algurithia  caa  easily  be  generalized  to  higher  space  dimen- 
sions. 

We  report  on  numerical  experiments  with  our  method.  We  solve  test  problems  with 
known  solutions,  and  we  compute  a  flow  evolving  from  an  initial  vortex  blob.  Our 
results  indicate  that  the  method  is  convergent  of  first  order. 

As  an  application  of  our  method,  we  present  a  fully  Lagrangian  variant  of  Peskin's 
algorithm  for  the  treatment  of  elastic  boundaries  immersed  in  the  fluid,  and  we  report 
on  numerical  results  obtained  with  this  algorithm. 


Acknowledgement!  .   /..  j 

I  am  deeply  grateful  to  Professor  Charles  S.  Pesldn  for  his  instructive,  enjoyable 
and  encouraging  way  of  working  with  me,  without  ever  becoming  impatient  at  my 
weaknesses. 

I  am  eqtially  grateful  to  Professor  Olof  B.  Widlimd  for  three  years  of  invaluable 

help  and  advice.  I  thank  him  for  aU  I  have  learned  from  mm,  for  all  the  time  and  con- 

.. .  "i^anos  tc3."isi  . . . 
sideration  which  he  has  devoted  to  my  r  ^ucation,  and  for  his  generous  kindness. 

Studying  at  the  Courant  Institute  has  been  a  great  and  rewarding  experience.  I  am 

thankful  for  all  I  have  been  taught,  for  the  friendly  and  generous  atmosphere  at  the 

"  '.  .1'.'. '  7.--  I : 
institute,  and  for  financial  support  from  the  Courant  Institute  and  the  Graduate  School 

of  Arts  and  Science  at  New  York  University. 

I  thank  the  mathematics  group  at  the  Lawrence  Berkeley  Laboratory  for  the 
interesting  and  pleasant  summer  which  I  spent  there  as  a  visitor  in  1984. 


1.  Introdactlon 

Most  of  the  methods  presently  used  in  computational  fluid  dynamics  use  the 
Eulerian  cather  than  the  Lagrangian  coordinate  system.  The  strength  of  Lagrangian 
methods  for  the  Navier-Stokes  equations  is  the  natural  and  simple  way  in  which  the  non- 
linear convc-tJ'M  ;^art  of  the  equadonj  ic  treated,  by  moving  the  grid  in  the  flow.  How- 
ever, this  roa/em-mt  of  the  grid  ca<j:ses  a^icial  problems.  The  discretization  of  the  dif- 
ferential operators  and  the  numeiic  fjlution  of  the  resulting  discrete  problems  are 
more  difficult  on  irregular  grids  )k^n  on  the  regular  grids  used  for  Eulerian  computa- 
tions. 

Thc^  difflcultifs  have  hampeied  the  development  of  fully  Lagrangian  methods  for 
flow  with  large  deform&tions  in  two  or  three  dimensions.  Methods  which  combine 
Lagrangian  and  Eulerian  fr^'ii  ■!;«  htfv^  be«n  studied  for  some  time.  Examples  are  the 
MAC  (Marker  and  Cell)  methrj-i  c*  io  i>j]ow  and  Welch  (1965),  in  which  Lagrangian 
markers  are  used  to  dr.ter.Aa^'i  the  •>OL-ili'jr  of  a  free  boundary,  but  not  for  the  solution 
of  the  Navier-Stokes  eqiiatioiffl?,  and  th-i  AT..P  (Arbitrary  Lagrangian-Eulerian)  method, 
in  which  rezoning  steps  fnevfnt  'i.e  mesh  from  changing  its  initial  topology;  src  Hlrt, 
Amsden  and  Cook  (1974),  Chan  (1975). 

There  has  been  recent  progress  in  overcoming  the  difficulties  inherent  in  the  use  of 
fully  Lagrangian  methods.  Discretizations  of  differential  operators  on  irregiilar  grids  can 
be  computed  with  increased  efficiency  by  using  new  computational  geometry  algorithms. 
The  solution  of  the  resulting  discrete  problems  has  become  easier  due  to  new  numerical 
techniques  such  as  multigrid  methods.  We  briefly  discuss  these  two  developments  in 
turn. 

In  order  to  discretize  a  differential  operator  on  a  general  set  S  of  points  in  the 
plane,  one  usually  needs  a  mesh  associated  with  S,  for  example  a  set  of  polygons  each 
containing  exactly  one  of  the  points  in  S,  or  a  triangulation  whose  vertices  are  the  points 
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in  S.  The  Voronoi  diagram  associated  with  S  can  serve  as  a  polygonal  mesh  or  as  a  tool 
for  the  construction  of  a  triangular  mesh.  The  Voronoi  polygon  associated  with  2C^5  is 
the  set  of  points  at  least  as  close  to  2C  as  to  any  other  point  in  S;  see  Voronoi  (1908). 
Recently  many  algorithms  for  the  construction  of  Voronoi  diagrams  have  been 
developed;  sec  Shamos  and  Hoey  (1975)  :for  a  very  fast  sequential  algorithm,  and 
Aggarwal,  Chazelle,  Guibas,  O'lXmlaing  and  Yap  (1985)  for  a  parallel  algorithm.  Pes- 
Idn  (1979)  introduced  an  algorithm  for  uptHting  Voronoi  diagrams  within  a  time  depen- 
dent calculation.  For  N  points,  0(N)  operations  are  required  per  time  step.  This  algo- 
rithm can  be  generalized  to  three  dimensions  with  few  changes.  We  note  that  Voronoi 
diagrams  have  been  used  in  Lagrangian  fluid  dynamics  codes  by  Augenbaum  (1982), 
Peskin  (1979),  Trease  (1981).  ,';iifir 

In  addition,  efHdent  iterative  reconnection.' algorithms  for  triangular  meshes  have 
been  introduced;  see  Crowley  (1971),  PHIkp' ''-.  C1981),  Fritts  and  Boris  (1979).  The  tri- 
an^ations  thus  generated  are  closelj^xelated  to  Vortmoi.  diagrams  and  can  in  fact  be 
xised  to  construct  Voronoi  diagrams.  The  generalization  of  these  algorithms  to  three 
dimensions  is  possible,  but  complicated;  see  Fritts  (1985). 

A  further  mesh  construction  algorithm  applicable  to  Lagrangian  fluid  dynamics  has 
recently  been  introduced  by  Boris  (1985).  Given  N=N^-N2  arbitrary  points  in  the  plane, 
this  algorithm  finds  an  ordering  2C(«i.«2),  l^ii^^i,  Is/j^^j  such  that 

This  indexing  is  called  an  MLG  (Monotonic  Logical  Grid).  It  is  easy  to  see  that  an  MLG 
always  exists  and  can  be  computed  in  0(NlogN)  operations.  Boris  describes  a  simple 
iterative  proceduie  by  which  it  can  be  updated  widiin  a  time  dependent  computation  in 
only  0(N)  operations  per  time  step.  The  algorithm  remains  unchanged  in  arbitrary 
dimensions  and  is  easily  vectorizable. 

We  turn  to  the  numerical  solution  of  equations  on  Lagrangian  grids.  Many  numer- 
ical  algorithms  for  incompressible  fluid  dynamics  require  the  solution  of  Poisson 


equations  for  the  pressure.  An  example  is  the  projection  method  introduced  by  Chorin 
(1968).  We  shall  consider  a  Lagrangian  version  of  the  projection  method  which  requires 
the  solution  of  a  Poisson  equation  for  the  pressure  and  a  Helmholtz  equation  for  each 
velocity  component.  We  use  a  multigrid  algorithm,  more  precisely  a  twogrid  algorithm, 
which  works  well  even  on  random  grids.  A  different  successful  multigrid  algonihm  for 
discrete  Poisson  equadons  on  irreguJar  grids  has  been  introduced  by  L5hner,  Morgan, 
Peraire  and  Zienkiewicz  (1985).  The  n^w  AMG  (Algebraic  Multigrid)  algorithms  appear 
to  be  still  another  possibility;  sec,  t.g.,  StQben  (1983).  They  have  apparently  not  yet 
been  used  in  Lagrangian  fluid  dynamics. 

Most  of  the  recent  work  on  Lagrangian  methods  concentrates  on  invisdd  flow.  In 
this  thesis,  however,  we  shall  study  a  Lagrangian  method  for  the  viscous,  incompressible 
Navier-Stokes  equations,  using  primitive  variables.  We  now  give  an  introduction  to  this 
method. 

In  1979,  Peskin  introduced  a  new  Lagrangian  method  for  the  incompressible 
Navier-Stokes  equations,  using  finite  difference  operators  constructed  with  the  aid  of 
Voronoi  diagrams.  The  method  requires  the  solution  of  a  linear  symmetric,  indefinite 
system  of  equations,  essentially  a  discretization  of  the  Stokes  equations,  in  every  time 
step.  No  efficient  solver  was  found  for  this  problem.  Only  few  numerical  tests  were 
therefore  performed  with  the  method,  all  using  very  small  grids. 

In  this  thesis,  a  fractional  step  method  which  is  a  modification  of  Pesldn's  original 
method  is  developed.  Only  discrete  Helmholtz  and  Poisson  problems  need  to  be  solved 
in  every  time  step.  The  much  improved  efficiency  due  to  the  two-level  iteration  has 
allowed  us  to  conduct  experiments  on  grids  of  substantially  larger  size.  The  present 
method  cannot  yet  compete  with  existing  Eulerian  methods  for  the  incompressible 
Navier-Stokes  equations.  We  believe,  however,  that  this  work  shows  that  the  method  is 
promising. 
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We  give  a  short  summary.  The  two-dimensional  incompressible  Navier-Stokes 
equations  have  the  form 

IC+U-^M-vAjt+Sp  =/  (1.1) 

Y-Jt  =  0  , 

where  ii=itis.,t)^R^,p=p(i.,t)(iR,S.iR^,t^R,  t>0.  We  impose  an  initio  condition 

M(i,0)=aoU)         -  •  ^^-  (1.2) 

and  the  periodicity  conditions  r      .  - 

it(i+i,0=itU,0.    pU+t,t)=pU,t),    t^Z^.  (1.3) 

In  section  2,  we  describe  a  way  of  adapting  Pesldn's  1979  algorithm  for  the  con- 
struction of  Voronoi  diagrams  to  the  case  of  a  periodic  domain.  Pesldn  also  introduced 

.;  i  '-it-  .  T 

discretizations  of  the  Laplace,  gradient  and  divergence  operators  on  irregular  grids  in  the 
plane,  based  on  the  Voronoi  diagram;  see  Pesldn  (1979).  We  present  the  derivation  and 
properties  of  these  operators  in  section  3i  In  section  4,  we  describe  the  two-level  itera- 
tion  and  give  results  of  numerical  experiinents  which  indicate  that  the  method  is  almost 
as  efficient  as  multigrid  algorithms  for  Helmholtz  equations  on  regular  grids.  In  section 
5,  we  consider  the  problem  of  projecting  a  vector  field  onto  its  divergence  free  com- 
ponent. We  describe  an  iterative  algorithm  whidi  is  applicable  for  irregiilar  grids  and 
which  quickly  finds  a  good  approximation  to  the  solution  of  the  continuous  problem. 
Ojnvergence  to  the  solution  of  the  discrete  projection  problem  is  provable  but  slow.  In 
section  6,  we  combine  these  tools  to  obtain  a  fractional  step  method  for  the  Navier- 
Stokes  equations.  Numerical  experiments  suggest  that  the  method  is  of  first  order  in  the 
sense  that  the  discretization  error  is  roughly  reduced  by  1/2  if  the  number  N  of  fluid 
markers  is  multiplied  by  4  and  the  time  step  is  reduced  by  1/2.  (In  two  space  dimen- 
sions, the  effective  meshwidth  h  is  0(A^~^),  therefore  multiplication  of  N  by  4  reduces 
h  by  2.)  The  work  per  time  step  is  0(NlogN)  and  could  be  reduced  to  0(N)  by  a  trivial 
modification.  The  method  is  implicit,  and  there  is  no  stability  condition  on  the  size  of  the 
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time  step.  In  section  7,  we  describe  an  algorithm  for  immersed  boimdaries  due  to  Pes- 
kin  (1977).  We  note  that  (1.3)  is  an  appropriate  condition  in  this  context.  The 
immersed  boundaries  are  represented  by  chains  of  Lagrangian  particles  connected  by 
springs.  For  the  fluid,  Pesldn  uses  an  Eulerian  method,  which  we  replace  by  our  Lagran- 
gian method.  We  present  some  preliminary  numerical  results.  In  section  8,  we  discuss 
difficulties  with  our  method  and  suggest  ways  of  overcoming  them. 
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2.  Construction  of  periodic  Voronoi  diagrams 

In  section  2.1,  we  outline  the  algorithm,  omitting  details  of  implementation  unless 
they  are  useful  for  a  general  description  of  the  method.  In  section  2.2,  we  give  algo- 
rithmic details  and  discuss  our  implementation  of  the  method. 

2.1  The  algorithm 

Let  5  be  a  set  of  points  in /?^.  For  2C€ 5,  we  define 

P(Z„S)  :=  U:€/?2  .  ||j._2^||  s  ii^.^ii  foy  jji  ^^sy  (2.1) 

P(Ili,S)  is  an  intersection  of  half  planes,  i.e.  a  convex  polygon.  It  is  called  the  Voronoi 
polygon  (Voronoi,  1908)  or  Dirichlet  region  of  X.  with  respect  to  S.  The  collection  of  all 
P(.X.,S),  2C€5,  is  called  the  Voronoi  diagram  associated  with  S  and  is  denoted  by  V(5). 
An  edge  in  V(5)  is  generated  by  two  points,  i.e.  it  belongs  to  two  polygons.  A  comer  is 
generated  by  three  or  more  points,  i.e.  it  belongs  to  three  or  more  polygons.  We  always 
assume  that  every  comer  is  generated  by  exactly  three  points.  A  degenerate  case  such  as 


is  resolved  into 


This  resolution  of  multiple  comers  into  sets  of  single  comers  linked  by  edges  of  length 
zero*  is  always  possible  but  never  unique.  This  fact  is  of  some  importance  for  our  con- 
struction algorithm  described  below. 

It  is  well  known  that  the  Voronoi  diagram  for  a  set  of  N  points  in  the  plane  can  be 
constructed  in  0(NlogN)  operations  and  that  this  operation  count  is  best  possible;  see 


Shamos  and  Hoey  (1975).  However,  if  additional  infc^aiation  on  5  is  available,  better 
operation  counts  can  be  achieved.  The  algorithm  introduced  by  Peskin  (1979)  makes  use 
of  the  assumption  that  the  Voronoi  diagram  has  already  been  constructed  for  a  set  5 
from  which  S  can  be  obtained  by  slightly  displacing  each  point,  and  that  certain  pj>ces  or 
information  about  this  construction  have;  been  saved.  Similarly,  the  algorithiu  useu  by 
Fritts  (1981)  can  br  viewed  as  an  algcrirhm  which  uses  V(5)  to  compute  V(5).  Experi- 
ments indicate  that  both  algorithms  construct  V(S)  in  linear  time.  If  5  is  a  random  set, 
there  may  be  algorithms  for  which  the  expected  value  of  operations  is  smaller  than 
0(>flogN).  This  is  particularly  plausible  if  some  information  is  given  about  the  probabil- 
ity distribution  of  S.  Numerical  experiments  supporting  this  statement  will  be  presented 
later  in  this  section. 

We  now  consider  tlie  case 

S  =  5v  •=  (Y^^t  •  1<;<^.  tez^ },  with  i^qo,l)2.  (2.2) 

In  our  fractional  step  method,  thr  2[/  will  tr  fluid  markers.  The  corresponding  Voronoi 
diagram  can  be  viewed  as  a  periodic  Voronoi  diagram  in  the  plane  or  as  a  Voronoi 
diagram  on  the  torus 

T:=Ryz^  (2.3) 

with  the  metric 

dax):=ttmi{  IlX-I+tll  :  k^Z^  }•  (2.4) 

We  describe  a  modification  of  the  algorithm  due  to  Peskin  (1979)  for  the  construction  of 

V(5v). 

We  store  information  about  the  comers  in  V(5v)  which  lie  in  [0;1)^. 

Proposition  2.1:   V(5v)  has  exactly  M=2N  comers  in  [0;1)^ 

Proof: 

One  could  use  the  fact  that  the  Euler  characteristic  of  the  torus  is  0.  Notice,  how- 
ever, that  the  Voronoi  polygons  on  the  tonis  are  not  necessarily  homeomorphic  to  disks. 
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This  is  shown  by  the  following  example,  in  which  the  two  polygons  are  homeomorphic 

to{2C€i?2:is|lx||s2}. 

-Voronoi  polygon  of  2Ci 


Therefore  an  additional  argument  is  needed  to  show  M=2N.  We  note  that  those 
polygons  on  the  tonis  which  are  not  homeomorphic  to  disks  also  cause  difficulties  in  our 
algorithm  for  the  construction  of  V(5);  sec  below. 

Our  proof  that  M=2N  is  based  on  the  fact  that  the  Euler  characteristic  of  the 
sphere  is  2.  Consider  m  by  m  copies  of  the  points  X.i,"'J^- : 


(m=3  here) 


Z' 


• 

%        m 
• 

• 

• 

■'A 

•       • 
• 

the  unit  square 

Consider  the  corresponding  Voronoi  diagram.  Intersect  it  with  [0;  m]^.  The  result  is  a 
polygonal  mesh.  We  extend  this  mesh  to  one  on  the  sphere  /J^ui*}  by  introducing  edges 
which  connect  the  four  comers  of  [0;  m]^  with  ».  This  introduces  a  degenerate  comer  at 
«  which  we  resolve  into  two  non-degenerate  comers  connected  by  an  edge  of  length  0. 
The  mesh  on  the  sphere  has  m^A^+4  faces.  Since  every  edge  joins  2  comers  and  every 
comer  belongs  to  3  edges,  the  nimiber  of  comers  is  2/3  times  the  number  of  edges. 
From  Euler's  formula,  we  now  conclude  that  the  total  nimiber  of  comers,  including 
those  at «,  is  2m2^+4. 

The  (m-2)  by  (m-2)  inner  squares  are  identical  up  to  translation  and  contain 
(m— 2)^A/  comers.  For  fixed  N,  the  number  of  comers  in  the  4m-4  boundary  squares  is 
0(m).  Therefore  we  have 


(m-2)2il/  +  0{m)  =  Inp-H  ■¥  2 


Letting  m  -  «,  we  obtain  M  =  2N.  ■ 

It  follows  from  this  that  there  are  exactly  3N  edges  (modulo  Z^)  in  V(5).  This  is  of 
interest  for  our  method  because  edges  in  V(5)  correspond  to  nonzero  couplings  in  the 
Hnite  difference  operators  presented  in  section  3. 

We  call  the  2N  comers  in  [0;1)^  the  stored  comers  and  call  the  points  2Ci,"-,  2C\-  the 
stored  points.  For  the  j-th  stored  comer,  we  store  its  cartesian  coordinates 
(XC(j),YC(J))  and  the  square  RAD2(J)  of  its  radius,  whidi  is  defined  as  its  distance 
from  its  generating  points.  We  introduce  the  notation  rcui(c)  for  the  radius  of  the  comer 
£.  The  cartesian  coordinates  of  a  comer  in  V(5)  can  be  written  in  the  form 

XC(IQ+KX,  YC(IQ  +  KY 
with  uniquely  determined  integer  indices  IC,KX,KY.  We  call  IC  the  basic  index  and 
KX,KY  the  shift  indices.  We  call  the  comer  with  basic  index  IC  and  shift  indices 
KX.KY  the  comer  (KX,KY,IC).  We  use  analogous  conventions  for  the  points  in  S,  i.e. 
the  fluid  markers  and  their  translations.   We  store  integer  indices 

IPT(j,nJ)  ,    »■  =  1,2,3  ,    n  =  0J,2  and 
ICR(i,nJ)  ,    »■  =  1,2,3  ,    n  =  0,l,2, 

which  have  the  following  meaning.  The  point 

(IPT(«,1J),IPT(/,2J),IPT(/,0J)) 
is  the  i-th  generating  point  of  the  j-th  stored  comer,  and  the  comer 

aCR(i.lJ),ICR(/,2J),ICR(i.0j)) 
is  the  i-th  nJghbor  comer  of  the  j-th  stored  comer. 

We  now  outline  the  algorithm  for  the  construction  of  V(Sv). 

Algorithm  2.1: 

Step  1:  Construct  V(Sj). 
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Step2:For  j=2,...,N: 

Step  2.1:  Construct  V(5y_iuU,})  from  V(5,_J  . 
Step  2.2:  Construct  VCS})  from  V(5y_iuU/}). 

Thus  in  step  2.1,  we  add  the  stored  point  i/  without  its  periodic  images  to  the  otherwise 
periodic  Voronoi  diagram.  In  step  2.2,  we  then  add  all  the  periodic  images  of  Xj- 

Step  1  is  trivial.  The  central  part  of  the  algorithm  is  a  procedure  for  the  construc- 
tion of  V(5u{l})  if  I^  5  and  V(5)  is  known,  i.e.  step  2.1.  This  procedure  was  described 
by  Pesldn  (1979)  and  Bowyer  (1981)  independently.  For  completeness,  we  present  it 
here. 

A  comer  c  in  V(5)  is  called  broken  by  I  if  it  is  not  a  comer  in  V(S\j{X})-  Y.  breaks 
£  if  and  only  if  Y.  Ues  in  the  open  disk  around  c  with  radius  rad{c)-  X.  always  breaks  at 
least  one  comer.  This  is,  in  general,  gxxaranteed  if  Y.  lies  in  the  convex  hull  of  S.  To  see 
this,  notice  that  the  convex  huU  of  S  is  exactly  the  union  of  the  Delaunay  triangles,  i.e. 
the  triangles  which  are  obtained  by  joining  all  pairs  of  points  in  S  with  adjacent  Voronoi 
polygons: 


Y.  lies  in  one  of  the  Delaunay  triangles.  Consider  the  circumscribing  drde  of  this 
Delaunay  triangle.  Its  center  is  a  comer  c  in  V(5).  Snce  Y  lies  in  the  interior  of  this  cir- 
cle, it  breaks  £. 

Proposition  2.2:  The  set  of  broken  comers  is  connected  in  V(5). 
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Proof: 


This  proof  is  due  to  D.  Goldfarb  (unpublished).   The  comers  of  a  fixed  polygon  in 
V(S)  which  are  broken  by  X.  form  a  connected  set: 

I 


comers  broken  by  Z 


Walking  along  the  boundary  of  PiI.,S{j{X}),  one  visits  all  polygons  which  have  comers 
broken  by  X..  Whenever  dP(Y.,S{j{X})  crosses  an  edge  of  V(S),  this  edge  has  exactly  one 
endpoint  which  is  broken  by  Y..  This  endpoint  is  a  common  comer  of  the  two  polygons 
bordered  by  the  edge.  Therefoie,  a  walk  through  all  broken  comers  is  constructed  by 
foUowing  aP(i:,5u{I}).  ■ 


dPa,Su{l}^ 


We  can  therefore  find  all  broken  comers  in  0(1)  operations,  once  we  know  one  of 
them.  We  can  find  a  first  broken  comer  in  the  following  way.  We  first  choose  some 
point  ^iS  close  to  I.  Good  ways  of  choosing  X.  will  be  discussed  below,  since  the  effi- 
ciency of  the  algorithm  depends  mainly  on  this  choice.  We  determine  a  comer  c.  of 
P(X.tS)-  For  this  purpose,  we  use  an  array  IPTCR  with  the  following  meaning.  The 
comer 

(IPTCR(l,t),IPTCR(2,*),IPTCR(0,*)) 
is  a  comer  of  the  k-th  polygon.   Notice  that  the  pointers  IPT  and  ICR  introduced  above 
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point  from  comers  to  points  and  comers  and  are  therefore  different  from  IFTCR.  Be- 
ginning with  c,  we  then  conduct  a  breadth  first  search  through  the  graph  of  comers  of 
V(S)  until  a  broken  comer  is  found. 

Once  we  know  all  comers  broken  by  X.,  we  can  find  all  comers  of  the  new  polygon 

P(X.,S[j{^))  by  »ising  the  facts  that  eadi  such  comer  lies  on  an  edge  between  a  broken 

and  an  unbroken  comer  in  V(5)  and  that  on  each  such  edge  there  is  a  new  comer.  It  is 

useful,  but  net  necessary,  to  store  the  neighbor  comers  of  a  comer  c  in  counterclockwise 

order  and  to  number  the  generating  points  such  that  the  edge  connecting  c  with  its  i-th 

neighbor  comer  lies  opposite  to  the  i-th  generating  point: 

2 


This  ordering  is  the  most  important  specifically  two-dimensional  feature  of  our  code.  It 
facilitates  the  finding  of  the  generating  points  of  a  new  comer.  If  the  new  comer  lies  on 
an  edge  between  a  broken  comer  Ci  and  an  unbroken  comer  £2  in  V(5),  and  if  cj  is  the 
j-th  neighbor  of  c^  (j  €  {1,2,3}),  then  the  generating  points  with  indices  j-fl  (mod  3)  and 
j-t-2  (mod  3)  of  Ci  are  also  generating  points  of  the  new  comer.  The  third  generating 
point  of  the  new  comer  is  the  newly  introduced  point  X.'- 


Notice,  however,  that  these  two  points  could  also  be  obtained  by  Intersecting  the  sets  of 
generating  points  of  a  and  a  in  V(S),  regardless  of  &e  order  in  which  neighbor  comers 
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and  generating  points  are  stored. 

It  remains  to  describe  step  2.2.  To  simplify  the  notation,  we  take  j=N  and  define 


Fv:=/'(2Cv-,5v-iU{Xv})  and 


(2.5) 


Ps.=P(Kf;^s). 


(2.6) 


It  is  clear  how  to  get  P^-  from  P^: 

Ps=Psn([Xs;i  -  0.5; Xv.:  +  0.5]  x  [^v^  -  0.5; X.v^  +  0.5]),  (2.7) 

where  Xf,-  ^  and  X^-^  are  the  coordinates  of  2^-. 

We  first  assiime  toat  the  decision  which  comers  of  Pv  to  cut  off  can  simply  be 
made  by  using  (2.7),  by  checking  the  coordinates.  It  then  seems  easy  to  find  th.r  radii  of 
the  new  comers  in  V(5v),  their  neighbor  comers  and  their  neighbor  points,  as  indicated 
in  the  following  figure: 


To  find  the  neighbors  of  the  new 
comers,  follow  " "  and 


•  •  •  •  •  :  edges  in  ViS^-^J;- 


■■  ^P.v; 


-  —  -  -  :  translations  of  edges  of  Pvl ' 


dP. 


Therefore  the  description  of  step  2.2  seems  to  be  completed.  However,  while  this 
procedure  works  correctly  in  most  cases,  it  occasionally  breaks  down  because  of  round- 
ing errors.   We  show  a  typical  situation  which  leads  to  a  breakdown: 
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Here  N=2,  2Ci  =  (0.0,0.0),  i:  =  (0.0,0.5).  Py  should  have  exactly  6  corners,  two  of 
which  are  identical  up  to  translation  with  two  others.  In  any  finitely  accurate  arithmetic, 
the  computed  number  of  comers  of  P^-  may  be  4,5,6,7  or  8.  This  statement  is  explained 
by  the  following  figtire,  which  shows  a  way  of  cutting  which  results  in  7  computed 
comers  of  Pt,-. 


^n 


The  situations  which  cause  difficulties  are  cases  in  which  a  polygon  is  adjacent  to  its 
own  translation.  This  is  most  likely  to  happen  at  early  stages  of  the  construction.  It  may 
however  even  happen  in  V(5v)  with  arbitrarily  large  N,  as  when  all  points  2Ci,"  JC^-  lie 
on  a  line  parallel  to  one  of  the  coordinate  axes.  Even  if  a  polygon  is  adjacent  to  its  own 
translation,  difficulties  such  as  the  one  illustrated  above  occur  only  with  a  probability 
which  tends  to  zero  with  the  machine  epsilon.  It  is  desirable,  though,  to  have  an  algo- 
rithm which  will  work,  for  any  given  set  of  points,  if  the  machine  epsilon  is  sufficiently 
small  but  positive. 

We  therefore  have  to  cut  /y  ^  ^  more  careful  way  to  obtain  /'v  Le^  £  be  a  comer 
of  Pff,  generated  by  £,  n,  Xff  Let  £  be  a  comer  of  Pv»  generated  by 
£+(1,0),  2j.+  (l,0),2C^-.  The  following  figures  show  two  such  situations.  For  simplicity, 
we  use  examples  in  which  the  points  lie  on  a  square  grid.  In  practice,  breakdowns  do 
occur  even  in  irregular  configurations,  even  though  with  small  probability. 


i    a=£+(i,o) 
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01+ (1,0) 


i  h  ,^:^^'0) 


^ 


^ 


-H-k     ll+(l,0)- 


(.)  (b) 

Case  (a)  is  the  one  discussed  above  as  an  example  for  a  breakdown.  A  situation  such  as 
case  (b)  occurs,  for  example,  when  constructing  the  Voronoi  diagram  associated  with  the 
set 

{2Ci=(0,0),2:2=(0,|),2C3=(y.j)} 
Here  the  point  being  introduced  is  Xa,-  =  X.3  =  iir,-r)- 

To  decide  whether  c  is  broken  by  2C^--(l,0)  in  situation  (a),  we  consider  the  Voro- 
noi diagram  V({i;  a;  Xv-;  Xv  -  (1,0)}): 


^-fi) 


«> 


If  we  draw  the  diagram  like  this,  we  conclude  that  c  is  not  broken  by  Xv  ~  (1,0),  since 
it  is  a  comer  in  the  four-point  Voronoi  diagram.  Qearly,  we  can  use  the  same  diagram 
to  decide  whether  c  is  broken  by  Xv+(1,0).  We  conclude  that  it  is.  If  we  had  drawn  the 
short  edge  of  length  0  like  this: 

31 


*^-  fi) 


X^• 


we  wotild  have  concluded  that  c  is  broken  whereas  c  is  not.  In  any  case,  we  make  oppo- 
site decisions  for  c  and  c,  which  is  correct  here.  We  consider  case  (b)  now.  In  this  case, 
it  is  not  correct  to  make  opposite  decisions  for  c  and  c.:  Neither  c.  nor  c.  are  broken  here. 
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Again,  we  can  read  this  off  from  the  Voronoi  diagram  V({t;  a,;  2C^•;  Xv  -  (1,0)}): 


^-6) 


£ 

0 


This  motivates  the  following  algorithm  for  cutting: 

Algorithm  2.2: 

For     eadi     comer     c.     of     Pf,-,     dedde     whether     or     not     it     lies     in 

[^.v.i  ~  0-5;-X^.v,i  +  0-5]  X  [Xs^  -  0.5;  Xv^  +  0.5],  using  the  following  procedure. 

Decide  whether  c  is  closer  to  2Qv  ^^^^^  ^  ^•■'"4o»  4o=±(1.0)  or  4©  =  ±(0,1)  by 
computing 

V(a;3i:2Qv;2Cv+4o}).  (2.8) 

where  Xv>£jl  are  the  generating  points  of  £.  If  £,  is  a  comer  in  this  diagram,  it  is 
at  least  as  dose  to  2Cy  as  to  2C^•  +  io>  otherwise  it  is  broken  by  <^-  +  Jq.  If  c  is  a 
comer  of  Py  generated  by  i^•  £  ~  io>ll  ~  Ao»  *^^  ^^  same  diagram  to  decide 
whether  or  not  c.  is  broken  by  2^-  -  io- 

In  our  code,  this  procedure  is  simpliHed.  The  main  simplification  is  a  test  whether 
cutting  is  necessary  at  all.  This  can  safely  be  done  by  checking  the  coordinates  of  the 
comers  of  Pv-  ^  ^^  ^^^  majority  of  cases,  no  cutting  is  necessary.  If  cutting  is 
reqxiired,  we  need  to  determine  which  of  the  four  unordered  triples  of  points  in 
{iai'2CN-'JC\-  +  4o}  generate  comers  in  the  diagram  (2,8).  This  is  easy,  since  it  is  simple 
to  dedde  for  any  given  triple  whether  it  forms  a  comer  in  (2.8).  There  are  two  or  three 
comers  in  (2.8),  as  shown  in  the  foUowing  figures. 
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Computiog  the  topological  structure  of  the  Voronoi  diagram  is  an  ill-conditioned 
problem  in  the  sense  that  the  solution,  the  topological  structure,  changes  discontinuously 
with  the  data,  the  given  points.  Therefore  the  difficulty  which  we  have  described  and 
overcome  should  not  be  called  numerical  instability.  We  rejected  the  first  version  of  the 
algorithm  not  because  it  computes  the  wrong  topological  structiire,  but  because  it  com- 
putes two  different  structures  at  tht  same  point  on  the  tonis.  The  improved  version  of 
the  algorithm  may  still,  due  to  roimding  errors,  compute  the  wrong  topological  structure: 


instead  of 


^ 


This  is  unavoidable,  but  it  fortunately  does  not  matter  in  our  application. 

We  remark  that  the  areas  of  the  Voronoi  polygons  depend  on  the  given  points  in  a 
continuous,  in  fact  in  a  differentiable  way;  see  section  3. 

We  finally  describe  two  ways  of  finding  a  good  starting  point  2C  for  the  seardi  for 
the  first  broken  comer.  If  the  points  in  J  are  known  to  be  roughly  uniformly  distributed 
in  [0;1]^,  we  can  divide  [0*4]^  into  n^  sqxiare  cells  of  equal  size  with  n*=A^  and  create 
pointers  froc  those  cells  to  the  points  in  S  which  they  contain.  A  linked  list  is  a  good 
data  structiire  for  this  purpose.  In  each  cell,  a  pointer  to  one  point  in  the  cell  is  stored. 
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Then  each  point  in  the  cell  points  to  another  in  the  cell,  with  the  last  point  containing  a 
stop  code.  We  introduce  the  points  cell  by  cell  in  the  order  indicated  in  the  following  fig- 
ure: 

n      n+1       ..       N-n+1 

..      n-t-2       ••       N-n+2 
2 
1        2n        ••  N 

When  introducing  2Cjt  >  ^  good  choice  of  X.  is  then  2L -i-  Using  this  method  on  random 
grids  in  [0;1]^,  we  obtain  the  following  CPU-times. 


Table  2.1.    CPU-times  per  point  in  msec  on  a  VAX  11/780,  20  uniformly  distributed 
random  grids  for  each  N. 


N 

worst  case 

best  case 

average 

400 

16 

13 

14 

800 

15 

13 

14 

1200 

15 

14 

14 

1600 

16 

14 

15 

2000 

16 

15 

15 

The  measured  CPU-times  support  our  assertion  that  the  CPU-time  per  point  is  bounded 
independently  of  N. 

In  our  fractional  step  method,  we  use  this  method  at  time  0.  At  later  times,  we  use 
the  following  method.  Assume  that  V(S)  has  been  constriicted  for  a  set 
S-^  +  k:j^{h---J^},  k^Z^},  where  the  distances  dQLj^)  are  small,  and  i,  ^  [0;1)2. 
Here  d  is  defined  as  in  (2.4).  Let  j{k^  be  the  index  of  a  point  in  S  whose  polygon  had  a 
broken  comer  when  ^  was  inserted  diiring  the  constnictioQ  of  V(5).  A  suitable  transla- 
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tion  of  the  point  ij^y.  is  then  close  to  2L  •   The  indices  ;(*o)  €  {I;-;  *o}  are  therefcre 

saved  while  constructing  V(S).  In  our  application,  5  is  the  set  of  fluid  markers  and  their 
periodic  images  at  a  given  time  step,  and  S  is  the  corresponding  set  at  the  preceding  time 
step. 

This  method  is  usually  slightly  more  efficient  than  the  first  one.  It  also  has  t!is 
advantage  of  not  reqxiihng  any  renumbering  of  the  points. 

We  conclude  this  section  with  scms  remarks  about  other  possible  ways  of  con- 
structing periodic  Voronoi  diagrams. 

It  is,  of  course,  possible  to  reduce  the  problem  to  nonperiodic  Voronoi  diagrams  in 
the  plane.  Consider 

i  :=  {2C/  +  t:  1  s  7  ^  ^,  iCZ^  |*J  s  1  and  pt^l  ^  1}.  (2.9) 

It  is  easy  to  see  that 

PiXjJ)  =  P(Kj^)  for  all ;.  (2.10) 

Since  S  contains  9N  points,  this  procedixre  is  quite  inefficient.  The  definition  of  S  can  be 
modified  to  reduce  the  number  of  auxiliary  points,  but  the  algorithm  then  becomes  com- 
plicated. It  becomes  impossible  to  search  for  the  first  broken  comer  using  the  second 
strategy,  since  the  set  of  auxiliary  points  becomes  dependent  on  &e  given  configuration 
and  therefore  time  dependent. 

It  may  also  be  possible  to  work  with  the  torus  (2.3)  and  the  metric  (2.4)  directly. 
The  point  insertion  algorithm,  however,  cannot  be  used  in  literally  the  same  way  as  in 
the  plane.  To  see  this,  consider  the  step  in  which  the  coordinates  of  a  new  comer  are 
computed  after  its  three  generating  points  2C,Z>Z  have  been  foimd.  In  the  plane,  this  is 
done  by  computing  the  midpoint  of  the  drde  through  i  J^,Z-  On  the  torus  T  with  metric 
d,  there  is  no  such  unique  midpoint,  as  is  shown  by  the  following  example. 
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• • 


Each  point  •"  has  the  same  distance  from  each  point  V> 

A  periodic  version  of  the  algorithm  used  by  Fritts  (1981)  appears  to  be  an  excellent 
alternative.  This  algorithm  constructs  the  Delaunay  triangulation  rather  than  the  Voronoi 
diagram  and  seems  more  efficient  than  our  method. 

However,  we  believe  that  an  approach  like  the  one  described  here  may  be  more 
generally  applicable  because  it  makes  use  of  the  metric  only  and  does  not  require  an 
inner  product  or  angles.  For  a  version  on  the  surface  of  a  sphere  see  A'.igenbaum  and 
Pesldn  (1985).  A  generalization  to  three-dimensional  periodic  Voronoi  diagrams  seems 
straightforward. 

2.2  Some  (tirther  algorithmic  details  and  implementation  of  the  method 

We  describe  a  Fortran  77  subroutine  DIAGRAM  which  implements  the  method 
presented  in  section  2.1.  The  development  of  this  subroutine  began  with  a  program  writ- 
ten by  C.  Peskin  for  the  construction  of  non-periodic  Voronoi  diagrams.  Each  of  the  fol- 
lowing sections  corresponds  to  one  of  the  larger  subroutines  called  by  DIAGRAM. 


2.2.1  Reordering  fluid  markera  at  time  zero 

We  use  integer  indices 

ICZ(»,i),  i=l,...^,  *=0,1,2. 
At  t(=time)=0,  these  indices  have  the  following  meaning.  When  inserting  the  i-th 

point,  the  search  for  a  broken  comer  begins  with  a  comer  of  the  polygon  of  the  point 


KCZ(OJ) 


ricz(i,i)^ 

ljCZ(2,i)J 
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See  sections  2.2.3,  2.2.5  for  the  meaning  and  use  of  ICZ  at  t>0.   At  t=0,  the  following 
algorithm  is  used  within  DIAGRAM  to  generate  ICZ. 

Algorithm  2.3: 

Step  1:  Choose  ns=V^.  Create  pointers  IBTOP(it,0€{l,...A},  (t,0€{l,...,np, 
and  JFTO?{i)(i{l,...J^},  je{l,...A},  with  the  following  meaning.  If 
IBTOP(it,/)  =  0,  then  there  is  no  fluid  marker  in  the  cell  [(Jk-l)A;  iWi]x[(/-l)A;  Ih], 
otherwise  IBTOP(it,/)  is  Ae  index  of  one  such  marker.  In  that  case,  the  list  of  all 
indices  of  markers  in  the  cell  is  i;,. ..,»,■  where  i,  =  IBTOP(Jk,/),  ij^.^=JPTOP(ij), 
and  j  =  r  is  the  first  integer  for  which  IPTOP(«y)  is  zero.  The  construction  of  these 
pointers  proceeds  marker  by  marker,  so  that  each  marker  is  associated  with  one 
and  only  one  cell,  even  if  it  lies  on  the  border  between  two  or  fo;vr  cells.  To 
reduce  the  required  amount  of  testing,  we  use  an  array  IN  of  size  nxn.  IN(k,l)  is 
the  number  of  markers  found  so  far  in  cell  (k,l). 

Step  2:  Sort  the  fluid  markers,  following  the  ordering  (k,l)  =  (l,l),...,(l,n), 
(2,n),...,(2,l),  (3,1),...   of  the  cells,  with  arbitrary  ordering  within  each  cell. 

2.2.2  Initialization  of  variables 

When  introducing  the  first  fluid  marker,  two  comers  with  the  same  coordinate  pair 
are  generated.  Information  about  these  comers  is  stored  on  IPT,  ICR,  XC,  YC,  RAD2, 
as  described  in  section  2.1. 

The  entries  IPTCR(0,1),IPTCR(1,1),IPTCR(2,1)  of  the  array  EPTCR  are  initial- 
ized. Generally,  at  a  stage  of  the  construction  at  which  iq  points,  iq^S,  have  been 
introduced,  IPTCR(m,j)  is  defined  for  m  =  0,l,2,   l^jtsj,,.  The  comer 

(IPTCR(l,0,IPTCR(2,i),IPTCR(0,j)) 
belongs  to  the  i-th  polygon  at  that  stage.  IPTCR(1,1),IPTCR(2,1)  may  be  =jtO. 
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A  logical  array  FLAGB  is  initialized  by  FLAGB(il,i2  j)=  .FALSE,  for  all 

(il,i2j)  e  {-i,o,ipx{i...Ar}. 

We  shall  describe  the  use  of  FLAGB  in  section  2.2.4. 

In  sections  2.2.3  •  2.2.9,  we  shall  give  a  detailed  description  of  the  insertion  of  all 
translations  of  the  k^-th  point  (2s)toSN)  into  the  diagram  associated  with  the  transla- 
tions of  the  first  itg- 1  points. 

2.2.3  Determination  of  the  comer  at  which  the  search  for  the  first  broken  comer  starts 

We  consider  the  case  t>0.  As  before,  we  use  the  notation  Xi,-- Jt\-  for  the  posi- 
tions of  the  points  during  the  previous  construction.  Let  iJll,•••,^l,^•C2^  be  such  that 
||(ik  +  U<i)-2Cill  is  small  (k=l,...,N).  Within  our  fractional  step  method,  the  pairs  y^  are 
determined  while  moving  the  fluid  markers.  If  Xi  moves  beyond  the  right  boundary  of 


[0;!]^,  for  example,  it  is  shifted  by 


0 


k        / 


ez^  such  that  the  shifted  point  lies  in  [O;!]^, 
and  jj'i,i:=  — 1.   UsuaUy,  v^  j^  will  also  eqiial  -1. 

ICZ  is  not  newly  generated  but  taken  as  generated  during  the  previous  call  of 
DIAGRAM.  It  has  the  following  meaning.  When  the  JtQ-th  point  Xt  was  introduced  dur- 
ing the  previous  time  step,  a  comer  of  the  polygon  of  the  point 

(ICZ(l,*o),ICZ(2,Jto),ICZ(0.*o)) 
was  broken. 

Let  Xj^  be  the  point  to  be  inserted  next.  The  basic  index  of  the  comer  at  which  the 
search  for  a  broken  comer  begins  is 

ICl=IPTCR(0,ICZ(0,*o)).  (2.11) 

The  shift  indices  are 

KXl=IPTCR(l,ICZ(0,*o))+ICZ(l,ito)-Hfii^4-Hicz(o^.i  (2-12) 
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and 

KYl=IPTCR(2,ICZ(0,ito))+ICZ(2,*o)+^Li^^-^i:c2(OJkJ^•  (2.13) 

From  (2. 12), (2. 13),  it  is  clear  that  the  pairs  y^  need  to  be  stored  only  temporaiij  They 
can  be  used  to  modify  ICZ(l,i),ICZ(2,Jt)  immediately  and  are  not  needed  again. 

Even  though  shift  indices  of  comers  broken  by  JC*  €[0;lj^  always  lie  in  {-1,0,1}, 

KXl  and  KYI  may  not  belong  to  {-1,0,1}.  The  comer  (KX1,KY1,IC1)  is  then  unlikely 
to  be  a  good  starting  point  for  the  search  for  a  broken  comer.  Lacking  any  better  possi- 
bility, we  redefine  IC1=1,KX1=0,KY1=0  if  CKX1,KY1)  <  {-1,0,1}%  thus  ensuring  that 
the  search  for  a  broken  comer  can  always  be  limited  to  comers  with  shift  indices 
between  -1  and  1.  Such  a  limitation  is  useful  because  the  number  of  storage  locations 
needed  for  the  array  FLAGB  (see  section  2.2.4)  is  then  18N,  which  is  the  number  of 
comers  in  nine  copies  of  the  unit  square. 

2.2.4  Search  for  ■  first  broken  comer 

Beginning  at  (KX1,KY1,IC1),  a  breadth  first  search  tlirough  the  graph  of  comers  is 
conducted  to  find  a  first  broken  comer.  The  array  ICR  contains  the  adjacency  structure 
of  this  graph.  We  state  the  algorithm. 

Algorithm  2.4: 

Test  whether  the  comer  (KX1,KY1,IC1)  is  broken.  If  so,  the  algorithm  terminates. 
Otherwise:  Set 

Kl=l,  LISTB(0,1)=IC1,  USTB(1,1)  =  KX1,  USTB(2,1)  =  KY1. 

USTB  will  be  the  list  of  tested  comers.  Set  FLAGB(KX1, KYI, IC1)  =  . TRUE.. 
Generally,  FlAGB(il,i2,j)  will  be  .TRUE,  if  and  only  if  comer  (il,i2,j)  has  been 
tested. 
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Repeat  the  following  steps: 

IC=LISrB(0,Kl),  KX=LISTB(1,K1),  KY=LISTB(2,K1).  For  each  of  the  neigh- 
bors-(KX2,KY2,IC2)  of  the  corner  (KX.KY.IC).  test  whether  KX2C{- 1,0,1}  and 
KY2g{- 1,0,1}  and  FLAGB(KX2,KY2,IC2)=. FALSE..  If  all  these  conditions  are 
satisfied,  test  whether  the  comer  (KX2,KY2,IC2)  is  broken.  If  so,  the  algorithm 
terminates,  otherwise  (KX2,KY2,IC2)  is  appended  ro  LISTS,  and 
FLAGB(KX2,KY2,IC2)  is  set  to  .TRUE..  After  all  neighbors  of  the  comer 
(KX,KY,IC)  have  been  tested  in  this  way,  Kl  is  increased  by  1. 

Eventually,  this  algorithm  will  dearly  find  a  broken  comer. 

2.2.5  Search  for  all  broken  comen 

Once  a  first  broken  comer  (KX1,KY1,IC1)  has  been  found,  the  indices  of  the  first 
generating  point  of  this  comer  are  stored  as  ICZ(0,Jto),  ICZ(l,*o),  ICZ(2,ito).  All 
entries  in  FLAGS  are  set  back  to  .FALSE.,  using  the  list  LISTS.  Then  an  analogous 
search  for  all  broken  comers  is  performed,  using  LISTS  to  store  broken  comers  and  set- 
ting FLAGB(KX,KY,IQ=.TRUE.  if  (KX,KY,IQ  is  broken. 

2.2.6  Compotation  of  comers  of  P^ 

The  comers  of  P^  are  computed  as  described  in  section  2.1.  The  following  infor- 
mation is  generated. 

NCO  =  number  of  comers  of  A  . 

XCO,YCX)  =  vectors  of  length  NCX)  containing  the  coordinates  of  those  comers. 
Here  and  in  all  vectors  of  length  NCO  described  below,  the  ordering  of  the  comers 
follows  the  boundary  of  the  new  polygon  in  counterclockwise  direction.  The  use  of 
diis  ordering  could  easily  be  avoided.  This  is  important  because  there  is  no  obvious 
analogous  ordering  in  three  dimensions. 
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RAD20  =  vector  of  length  NCO  containing  the  squares  of  their  radii. 

rrO  =  integer  array.  For  y=l,---,NCO,  there  is  an  edge  from  a  comer  broken  by 

2L  to  the  unbroken  comer  (ITO(l, j).n"0(2,i), 110(0, j))- 

J) 

FLAGBO  =  boolean  vector  of  length  2N;  FLAGBO(j)=.TRUE.  if  and  only  if  the 
j-th  stored  comer  is  broken  by  some  translation  o(  X*  ,  i-t.   if  and  only  if  there  is  a 

pair  (il,i2)  such  that  FLAGB(iLi2,j)=.TRUE.   (ji{l,---,2kQ  -  2}  ). 

LISTBO  =  integer  vector  of  length  NBROKEN,  giving  the  basic  indices  of  all  bro- 
ken stored  comers. 

LISTP20,LISTP30  =  integer  arrays.  For  ;€{1,--,NC0},  the  generating  points  of 
the  j-th  comer  of  the  newly  constr\icted  polygon  P^  are,  in  counterclockwise  order, 

]^  and  the  points 


"0 


(USTP20(lJ),US'rP20(2J),USTP20(0j)), 
(USTP30(1J),USTP30(2J),USTP30(0J)). 

LISTLO  =  vector  of  length  NCO.  For  ;€{1,---,NC0},  the  LISTLO0>th  neighbor  of 
comer  (170(1, j),ITO(2,j),ITO(0,j))  is  broken  by  2Ci  ,  and  the  j-th  comer  of  P^  lies 

on  the  corresponding  edge. 

The  arrays  ICR,  IPT,  XC,  YC,  RAD2  which  describe  the  periodic  diagram  are  not 
yet  changed. 

2.2.7  Compotatlon  of  comers  of  P,^ 

First  we  test  whether  aU  (XC0(j),YCX)(j')),  j=l,...,NCO,  lie  in  the  unit  square  cen- 
tered at  Xi .  If  so,  no  comers  of  Pl  need  to  be  cut  off  due  to  interaction  with  periodic 

images;  see  section  2.1.  Otherwise  we  determine  the  comers  to  be  cut  off,  using  the  fol- 
lowing algorithm. 
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Algorithm  2.2b: 

Initialize  the  entries  of  an  integer  vector  IBR  of  length  NCO  by  IBR(j)=0  for  all  j. 
IBR0')=1  (2,3,4)  will  mean  that  the  j-th  new  comer  piC0(j),YCO(j))  is  broken  by 
the  northern  (western,  southern,  eastern)  translation  of  2L  •  It  is  easy  to  see  that 

no  other  translation  of  Xi   can  break  a  comer  of  P^ .  Initialize  the  entries  of  a 

boolean  vector  DONE  of  length  NCO  with  values  .FALSE.. 

For  j^{l,---^},  test  whether  there  is  a  Jq^{1,---J^}  such  that  the  set  consisting  of 
the  three  generating  points  of  the  comer  (XCO(j),YCOO'))  and  ^  -  fi]  is  equal,  up 
to  translation,  to  the  set  consisting  of  the  three  generating  points  of  the  comer 
(XCOOo) , YCOOo))  and  X*  +  (J) .  See  section  2. 1  for  illustrating  figures. 

If  such  a  JQ  does  not  exist,  test  whether  XC0<2Ci,i-y.  If  so,  set  IBR(j)=2, 
DONE0>.TRUE.. 

If  such  a  Jq  does  exist,  compute  the  corresponding  four-point  Voronoi  diagram  (see 
section  2.1)  and  set  DONEO)  =  DONE(;o)=.TRUE..  Set  IBR0')  =  2  if 
(XCOG),YCOG))  is  found  to  be  broken  by  2Ct,-(o),  ^^(Jo)=4  »f 
(XCOOo), YCOOo))  »  found  to  be  broken  by  Kt  +  (A). 

Finally,  test  for  all  ;€{l,-gv}  with  DONE0)=.  FALSE,  whether 
XC00)>2Ci,.i+  j;  if  50,  set  IBR0)=4. 

Set  all  entries  in  DONE  back  to  .FALSE,  and  search  for  comers  (XCOOO.YCOO)) 
broken  by  the  northern  and  southern  translations  of  2Ct  in  exactly  the  same  way. 

If  there  are  comers  with  JBR(J)=2,  there  are  comers  with  IBR(j)=4.  The  same  is 
true  for  comers  with  IBR0')=1,3;  see  the  following  figure. 
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IBR=1 


identified  with  each  other 


HS;^ 


identified  with  each  other 


mR=3 


Because  of  rounding  errors,  this  rule  might  be  violated.  Since  the  program  relies  on  it 
later,  it  sets  JBR(j)  from  i=l  (2,3,4)  back  to  0  if  there  are  no  indices  j  with  ffiR(j)  =  3 
(4,1,2).  This  does  not  cause  significant  errors  in  the  computed  shapes  of  the  polygons. 

After  having  defined  IBR,  it  is  easy  to  create  the  list  of  comers  of  the  polygon  Pl 

of  2Cji  in  the  new  periodic  diagram.  We  update 


"9 


NC0,XC0,YC0,LISTP20,LISTP30,RAD20,IBR,ITO,USTL0.  (2.14) 

NC0,XC0,YCX),LISTP20,LISIT'30,RAD20  are  updated  in  the  obvious  way.  IBR0')  =  O 
now  means  that  no  translation  of  X,..,  other  than  2Ct   itself,  is  among  the  generating 

points  of  the  comer  (XCOOO.YCOCj)).  roRO')  =  l  (2,3,4)  means  Lhat 

^^G)  (i.-(V).i.-(-°,).i/fi)) 

is  one  of  the  generating  points  of  (XCO(j),YCO(j)).  rrO(- J)  does  aot  have  any  mean- 
ing any  more  unless  IBR(j)  =  0,  nor  does  LISTLO  unless  IBR(j)  =  0  and 
FLAGB0(ITO(0J))=  .FALSE.. 

Comers  with  IBR=1  (4)  can  be  identified  with  comers  with  IBR=3  (2),  as 
explained  above.  Each  such  pair  represents  only  one  comer  in  the  new  periodic 
diagram.   Therefore  NCO  is  not  yet  the  number  of  comers  of  Xi  which  will  be  stored. 

To  prepare  the  insertion  of  the  new  comers  into  the  new  periodic  diagram,  indices 
IDENT(0,j),  IDENT(l,j),  IDENT(2,j),  j=l,...,NCO,  are  computed.  If  IBR0)^{1,4}, 
then  IDENT(0,j)  =  j,  IDENT(l,j)  =  0,  IDENT(2,j)  =  0.  If  IBR0>1,  then  IDENT(l.i)  =  0, 
IDENT(2,j)  =  l,   and  IDENT(0,j)   is  defined  as  indicated  in  the  previous  figure.   If 
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IBR04,  then  IDENT(l,j)=l,  IDENT(2,j)=0,  and  IDE^^^(0,j)  is  defined  analogous- 
ly- 

2.2.8  Compotation  of  edges  between  comers  of  Pl  and  other  comers 

There  are  two  cases  to  be  distingmshed  here.  If  IBR(j)  =  0  and  if 
FLAGB0(rrO(0,j))  is  .FALSE.,  then  the  comer  (ITO(l,j),rTO(2,j),ITO(0,j))  is  the 
neighbor  comer  of  the  j-th  comer  of  Pl  in  the  new  periodic  diagram.  If  IBR(j)^0  or 

FLAGB0(rrO(0,j))=.TRUE.,  the  neighbor  comer  of  the  j-th  comer  of  P^  in  the  new 

diagram  is  a  translation  of  one  of  the  comers  of  P^ ,  say  the  1-th  one.  This  translated 

comer  is  uniquely  characterized  by  the  fact  that  those  of  its  generating  points  given  by 
LISTP20,USTP30  are  the  points  whidi  are  also  stored  as  USTP20(-,j),USTP30(-,j).  It 
is  found  by  searching  through  all  comers  of  P^ .  We  perform  this  search  only  if 

IBR0")=2  or  IBR0')=3,  since  comers  with  IBR=1  or  IBR=4  will  be  identified  with 
comers  with  rBR=2  or  IBR=3. 

2.2.9  Insertion  of  the  new  polygons  into  the  periodic  diagram 

At  this  stage  of  the  algorithm,  LISTBO  contains  a  list  of  all  basic  indices  of  the  bro- 
ken comers.  Theses  indices  serve  as  indices  for  the  newly  created  comers.  Exactly  two 
additional  basic  indices  are  needed,  namely  2*0-1  and  2*o-  An  integer  vector  HISTBO 
of  length  NCO  is  aeated.  It  assigns  to  each  index  j^{l,...,NCX)}  in  the  circular  chain  of 
new  comers  a  basic  comer  index  in  the  new  periodic  diagram,  i.e.  one  of  the  indices  on 
USTBO  or  2*0-1,2*0.  Since  XCOOO.YCOO)  are  not  necessarily  ^[0;1],  XC(ILISTBOG)), 
YC(ILISTBOO'))  are  not  necessarily  the  same  as  XCOO.YCOO').  The  integer  differences 
XC(ILISTBOG))  -  XCOO),  YCCELISTBOO))  -  YCOQ)  are  stored  as  integer  arrays 
ISHX,ISHY.  Using  lUSTBO,  USTP20,  USTP30,  ICRO,  ISHX,  ISHY,  RAD20,  it  is 
then  easy  to  update  IPT,  ICR,  IPTCR,  and  RAD2,  completing  the  insertion  of  the  *o-th 
point  and  its  translations  into  the  periodic  diagram. 


-29- 

3.  Finite  difference  operators  on  Irregnlar  grids 

We  review  the  defmitions  of  discrete  Laplace,  divergence  and  gradient  operators 
due  to  Peskin  (1979)  and  discuss  their  properties.  Even  though  we  use  a  terminology 
appropriate  for  2  dimensions,  the  definitions  and  statements  in  this  section  are  dimen- 
sion independent. 

We  introduce  the  foUowing  notations.  For  2C€5v,  let  V[i]  be  the  area  of  the  Voro- 
noi  polygon  P(2C,5v)  =:P[K]-  For  If^v,  1*X.,  let  A[X.,I.]  be  the  length  of  the  common 
edge  of  P(.X.,S^)  and  P(I,5v).  A[K,r]  =  0  if  /'(2C,5v)  and  P(I,5v)  are  not  adjacent  to 
each  other.  Our  discrete  operators  couple  two  distinct  points  2C,1I  to  each  other  only  if 
they  are  neighbors,  i.e.  if  A[2i,]C]*0.  In  the  following,  we  let  Nb[ii]  be  the  set  of  neigh- 
bors of  2C.  By  definition,  2C^  Nb[X]. 


3.1  Discrete  Laplace  operator 

We  define  a  discrete  Laplace  operator  L  by 

V[2C](i:<l))(2C)  :=   2  A[2C.i:]-^^^^^^f^    for2C65v.  (3.1) 

This  definition  is  motivated  by  the  formula 

/A«^(i)4I  =  f^U)ds  (3.2) 

for  smooth  functions  <}>.  The  sum  in  (3.1)  is  formally  infinite  but  contains  only  fmitely 
many  non-zero  terms.  Similar  discretizations  of  the  Laplace  operator  have  been  used  in 
numerical  computations  for  a  long  time;  see,  for  example,  MacNeal  (1953).  The  follow- 
ing proposition  was  pointed  out  to  us  by  B.  Merder  (unpublished  communication). 

Proposition  3.1:  The  operator  in  (3.1),  seen  as  a  matrix,  is  the  stiffness  matrix  obtained 
with  piecewise  linear  finite  elements  on  the  Delaunay  triangtilation. 
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Proof: 

Both  matrices  have  zero  row  sums.  Therefore  we  need  to  consider  the  non-diagonal 
entries  .only.  The  matrices  have  the  same  non-zero  structiire  and  are  both  symmetric.  It 
therefore  makes  sense  to  talk  about  the  coupling  between  2C  and  2!  (  2C,21^Sv  i'^Z). 
We  compare  the  couplings  in  the  matrix  in  (3.1)  with  those  in  the  finite  element  stiffness 
matrix. 

Without  loss  of  generality,  we  assimie  that  2C  =  (0,0)  and  X.  =  (0,1).  Let  A,fl  be 
the  points  in  5v  such  that  AJLX  and  A  JC>J^  are  triangles  in  the  Delaunay  triangiilation: 


--^- 


I 
2C 


B. 
(solid:  Voronoi  diagram;  dotted:  Delaunay  triangles) 


Let  C  be  the  center  of  the  circle  through  AJC.I,  and  let  12  be  die  center  of  the  drde 
through  B.,X.,I.-  C  and  Q.  are  comers  in  V(Sf,)  •  The  claim  is  then  that  the  coupling 
between  X.  and  X.  obtained  with  piecewise  linear  finite  elements  is  [  C-fi  J . 

To  show  this,  we  first  compute  this  coupling  in  terms  of  the  coordinates  of  A  and 
fi.  The  relevant  finite  element  basis  functions  on  the  left  triangle  are 

<t>i(i)   =    ^^x,-X2+l,  (3.3) 


Ai 


the  basis  fimction  associated  with  2C,  and 


«|»2U)  =   -T-'i+'2.  (3-4) 

the  basis  function  associated  with  X.-  We  multiply  the  gradients  of  these  two  functions, 
integrate  the  product  over  the  left  triangle  and  multiply  by  -1.  The  result  is 


2Ai  2  ■ 


(3.5) 
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The  contribution  from  the  right  triangle  is  then,  by  symmetry, 

Bl-B,        B, 


(3.6) 


2^1  2  ■ 

The  difference  in  sign  between  (3.5)  and  (3.6)  is  due  to  the  fact  that  the  left  triangle  has 
the  area  -A/1,  while  the  right  triangle  has  the  area  +B/2.  We  have  then  shown  that 
die  finite  element  coupling  between  H  and  X.  i^ 


Al     ■"   A-l  At  0^'~Bn  By 

— --ir—^ -•¥—  .  (3  7) 

2Ai  2         2fli        2  ^      ' 


We  next  compute  |l£-I2| 


1  A-,(A-,-\) 

C  =  (y^i  +  -^^xr^'0-5)  (3-8) 

D.  =   (yfl:  +     -^^      .0-5)  (3.9) 

Therefore  the  square  of  (3.7)  rquals  the  square  of  |lC-d|. 

To  conclude  the  proof  of  our  assertion,  it  remains  to  be  shown  that  (3.7)  is  posi- 
tive. We  notice  that  the  foregoing  computations  do  not  make  any  use  of  the  fact  that  we 
are  using  Delaimay  triangulations  and  Voronoi  diagrams.  However,  all  off -diagonal  ele- 
ments of  the  finite  element  Laplacian  obtained  with  a  triangxilation  t  are  positive  if  and 
only  if  T  is  a  Delaimay  triangulation. 

We  outline  the  proof  of  this  well-known  result. 

A  triangulation  t  is  called  locally  equiangular  if  it  has  the  following  property.  If  T^ 
and  T^  are  adjacent  triangles  in  r  whose  union  Q  is  a  convex  quadrilateral,  the  sum  of 
the  angles  in  Q  which  are  cut  by  the  common  edge  of  7,  and  7;  is  larger  than  180". 


but  not 
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It  is  obvious  that  Delaunay  triangtilations  are  locally  equiangular.  Sibson  (1978) 
showed  that  the  converse  is  also  true. 

A  simple  calciilation  shows  that  all  off-diagonal  elements  of  the  finite  element 
Lapladan  obtained  with  a  triangulation  t  are  positive  if  and  only  if  t  is  locally  equiangu- 
lar; see  Fritts  (1981).  We  note  that  this  condition  is  also  equivalent  to  the  diagonal  dom- 
inance of  the  finite  element  stiffness  matrix,  since  this  matrix  has  the  zero  row  sum  pro- 
perty. 

This  concludes  the  proof  of  proposition  3.1.  ■ 

In  spite  of  proposition  3.1,  a  discretization  of  the  Poisson  eqxiation  based  on  (3.1) 
is  not  a  common  finite  element  discretization.  The  difference  lies  in  the  right-hand 
sides.  Let  /  be  a  given  continuous  right-hand  side.    In  the  finite  element  method,  the 

discrete  right-hand  side  consists  of  integrals  of  the  form  //(i)^(i)^,  where  /  is  an 

a 

approximation  of  /  and  ^  is  a  basis  function  of  the  finite  element  space.  When  using 
piecewise  linear  elements,  the  conmion  choice  of  /  is  the  piecewise  linear  interpolant  of 
/.  The  value  in  2C  of  the  discrete  right-hand  side  is  then  a  weighted  sum  of  values  of  /  in 
neighbors  of  2C-  The  sum  M[K.]  of  the  weights  (elements  of  the  mass  matrix)  is  one  third 
times  the  sum  of  the  areas  of  the  triangles  to  which  X.  belongs.  Using  (3.1),  one  is  lead 
to  using  V[X.]f(.X.)  as  the  value  of  the  discrete  right-hand  side  in  X..  V[K.]  and  A/[2C]  differ 
from  each  other  in  general,  but  they  are,  of  course,  equal  on  the  average.  M[X]  is  the 
area  of  a  cell  suggested  by  Dukowicz  (1981)  and  others,  which  is  obtained  by  joining  the 
centroids  of  the  triangles  with  the  midpoints  of  their  sides.  We  draw  the  two  different 
cells  in  one  figure: 


Voronoi  cell 


median— based  cell         '^A^' 

(The  median- based  cell  need  not  be  convex.) 
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We  now  prove  that  the  operator  L  is  weakly  consistent  to  first  order  with  the  con- 
tinuous Lapladan.  This  proof,  formulated  as  a  proof  of  weak  consistency  of  the  discrete 
divergence  operator  D  defined  in  section  3.2,  is  due  to  C.  Peskin  and  is  contained  in 
Bfirgers  and  Peskin  (1985).  We  define 


h  :=  max  diam(/»[i]) 


(3.10) 


Proposition  3.2:  If 


then 


(3.11) 


^H:i^)L6(x,)v[Xt]  =  f  Mi.)^Hi.)di  +  0(h) 

*  [0;!)' 

for  arbitrary  smooth,  periodic  scalar  functions  4>  and  v(». 


(3.12) 


Proof: 


We  approximate  the  difference  to  be  estimated  by 


(3.13) 


Since  <|>  is  uniformly  continuous,  the  error  which  we  have  introduced  is  only  0(A).  We 
next  apply  Green's  formula  (3.2)  and  its  discrete  counterpart  (3.1)  and  obtain 


2  2  •Kic*) 


A[^X] 


pwnpm 


da 


(i)ds 


(3.14) 


Since  the  expression  in  the  bracket  is  of  order  O(A-),  and  because  of  (3.11),  we  scsm  to 
obtain  an  upper  bound  of  size  0(1)  now.  However,  observe  that  the  roles  of  2^*  and  I  in 
(3.14)  may  be  reversed,  therefore  (3.14)  is  equal  to 

1 


-  2  2  ^U) 


^[i*,!] 


Hr)  -  Hx*) 

lll-2C.ll 


fwnpm 


3<t) 
da 


(i.)ds 


(3.15) 
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Here  a  always  denotes  the  normal  on  /*[2Ci]n''tt]  which  is  exterior  with  respect  to 
P[Xi]-  Since  (3.14)  and  (3.15)  are  equal,  they  are  also  both  equal  to  the  average  of  the 
two  expressions, 


T  2  2  (HX,)-Hr)) 

which  is  0(h).  ■ 


.    (3.16) 


Assumption  (3.11)  is  not  always  satisfied.  As  an  example,  consider  a  case  in  which 
aU  2Ct  lie  oo  a  vertical  or  horizontal  straight  line.  Then  A  =  0(1),  independently  of  N. 
(3.11)  is  a  non-degeneracy  assumption  on  the  polygons. 

We  give  a  second,  shorter  proof  of  (3.12).  Because  of  proposition  3.1,  (3.12)  is 
equivalent  to 

ai^„^,)  =  fl(<t.,4»)  +  0(A),  (3.17) 

where  a(-,-)  •*  ^^  bilinear  form 

fl(4»,«l'):=    J  l.Hs.)ZHi)di.  (3.18) 

associated  with  the  Laplace  operator  and  ^i,,^/,  are  the  piecewise  linear  interpolants  of 
(J>,»|;  with  respect  to  the  Delaunay  triangles.  (3.17)  is  well  known  from  the  theory  of 
interpolation  in  Sobolev  spaces  needed  in  finite  element  convergence  theory;  see  Ciarlet 
(1978).  The  non-degeneracy  assumption  on  the  triangulation  which  is  needed  for  this 
argtmient  is,  however,  neither  necessary  nor  sufficient  for  (3.11)  to  hold. 

We  gave  the  first  argument  for  two  reasons.  First,  it  is  independent  of  the  connec- 
tion with  the  finite  element  method  and  therefore  applicable  even  if  the  control  areas  are 
not  Voronoi  polygons  but  areas  defined  in  some  different  way.  Second,  it  is  applicable 
to  discretizations  of  differential  operators  other  than  the  Laplace  operator;  see  below  for 
the  case  of  the  divergence  operator. 
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Propositlon  3.3:  L  is  not  pointwise  consistent  with  the  Laplace  operator. 
Proof: 

We  give  a  counterexample.  Consider  the  following  Voronoi  diagram. 


(A*05 


We  use  Taylor's  theorem  to  expand  the  expression 

— li n —  V  17^  +  ~ii li —  V  Tz^  +  4>U2)-<^(io) 

llii-ioll  16  lll3-ioil  16 

around  xn.    The  coefficient  before  the  mixed  derivative  - — 2—  in  this  expansion  is 

—  — — A^,  which  proves  the  inconsistency  of  L.  Nothing  is  spedal  about  this  configura- 
tion. Almost  any  configuration  without  symmetry  properties  could  serve  as  a  counterex- 
ample. ■ 

L  is  pointwise  consistent  of  order  0  in  the  following  sense. 

Proposition  3.4:  For  a  fixed  smooth  periodic  function  4>,  L(t>  is  bounded  independently 
of  the  fluid  marker  configuration. 

Proof: 


V[2C] 


r»jc 
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for  some  constant  C  dqjcnding  on  (t>  only.  Now  notice  that 

2A[2C.I]lli:-2C||  =  4V[2C]. 


area  o?  this  triangl^i: 

This  concludes  the  proof  of  proposition  3.4.  ■ 

For  later  reference,  we  define  the  discrete  L^-product  by 


for  scalar  ftinctions  and 


i^,^)  :=  2<l>(2C*)4'(2Ck)V[it]  (3.19) 

k 


(M,3i)  :=  2ii(X*)li(2Ci)V[ii].  (3.20) 

k 


for  vector  fields. 


Also  for  later  use,  we  state  c\nd  prove  the  following  simple  algebraic  properties  of 
L. 

Proposition  3.5:   L  is  symmetric  and  negative  semidefinite  in  the  inner  product  (3.19). 
The  kernel  of  L  is  the  set  of  constants. 

Proof: 

(L^,^)  =  22  M^X]^^^l~ 'l^^UiiQ.  (3.21) 

Reversing  the  roles  of  ^  and  i.,  we  obtain 

(L^,^)  =  -2  2  A[^X]  ^^^l  ~  l^^^  ^ai  (3.22) 

kim^  ILL  -  AtW 


37 


Taking  the  average  of  (3.21)  and  (3.22),  we  obtain 


1  MZ)  -  «t)(2Ci) 

(L^A)  =  -y2  2  A[i„x.r'^^_2^\Hr)  -  hk,)),         (3.23) 


from  which  the  assertions  follow 


3.2  Discrete  divergence  operator 

We  turn  to  the  definition  of  the  discrete  divergence  operator, 

V[^p^U)  :=  2u(I)~^  for  2C€5v  .  (3.24) 

Here  — -^^  is  the  gradient  of  V[2C]  with  respect  to  H,  keeping  all  other  points  fixed;  see 

proposition  3.7.    We  also  use  the  notation  D  =:Di  to  distinguish  D  from  a   different 
discretization  Z?2  of  ^"  introduced  below. 

If  (2Ci,"- Jtv)  =  (2Ci(0."»^(0)  ^c  points  moving  at  velocity  u,  then 

j^VUjiO]  =  V[Xjit)]Diiaj(t))  (3.25) 

for  all  j,  in  particular: 

ProposidoD  3.6:    If  X.i(t),---Jijf(t)  are  points  in  space  which  move  in  a  discretely 
divergence-free  velocity  field,  the  Voronoi  polygons  maintain  their  areas  exactly. 

However,  the  Voronoi  polygons  do  change  their  areas,  usually  even  quite  drasti- 
cally, if  the  points  are  moved  in  a  continuously  divergence-free  field.  As  pointed  out  by 
J.  Dukowicz  (personal  communication),  this  is  even  true  for  large  numbers  of  fluid 
markers  and  small  diameters  of  the  Voronoi  polygons.  We  give  an  example  which 
proves  the  correctness  of  this  observation.  Consider  a  flow  of  the  form 

IL  =  M.  (3.26) 

with  a  constant  matrix  A.  Let  A  have  trace  zero,  so  that  (3.26)  describes  an  incompressi- 
ble flow.  One  could  think  of  (3.26)  as  a  local  approximation  to  a  real  flow  with  a 
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stagnation  point.  Let  a  set  of  fluid  markers  move  in  this  flow  and  consider  the 
corresponding  Voronoi  diagrams.  The  Voronoi  polygons  do  not  necessarily  maintain 
their  areas.   As  a  specific  example,  consider  a  flow  where 


A   = 


(0  U 
0  0 


(3.27) 


If  one  of  the  markers  is  at  the  origin  i=Q  at  time  0,  it  does  not  move.  Its  Voronoi 
polygon  at  later  times  depends  on  markers  which  are  far  away  from  it  at  time  0  and  will 
therefore  usually  not  maintain  its  area.  If  the  coordinates  of  all  fluid  markers  at  time  0 
are  multiplied  by  c,  then  the  resulting  Voronoi  diagrams  are  just  scaled  by  the  same  fac- 
tor e. 

In  summary,  this  example  shows  that  Voronoi  polygons  of  a  set  of  points  moving 
in  a  continuously  divergence  free  flow  may  change  their  areas  over  a  fixed  time  by  a  fac- 
tor which  does  not  converge  to  1  when  the  set  is  made  finer  in  the  sense  that  h,  as  in 
(3.10),  tends  to  zero.  Furthermore,  the  example  shows  that  the  same  statement  is  true 
for  all  control  areas  which  are  just  scaled  by  e  when  the  given  set  of  points  is  scaled 
by  €. 

An  explicit  formula  for  — ^^  was  derived  by  Peskin  (1979).  For  completeness, 
we  present  a  (slightly  different)  derivation  here. 

Proposition  3.7:  The  mapping  from  {(2Ci,-JCv)€(0;l)^'^'  :  i^*K.j  (oi  i*j}  into  R^' 
which  maps  the  cartesian  coordinates  of  2Ci,"  .Xv-  onto  the  areas  V[2Cj,- ••,V[2C\]  is  con- 
tinuously differentiable.  The  partial  derivatives  are  given  by 

ax.  Ill-ill    ^^'^^'  ^^^ 

11-^^^ II  is  not  bounded  for  r  ..  2C. 
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Proof: 

We  distinguish  the  following  three  cases  for  a  given  Z: 

i)2C  =  I 

n)X.^Nb[][](soX.*r) 
iii)2C^Iand2C^A^*[i:]. 
We  first  observe 

^^^  =  0  in  case  (iii)  .  (3.29) 

Case  (i)  may  be  reduced  to  case  (ii)  using  (3.29)  and  one  of  the  two  formulae 

2-^1^  =  0.  (3.30) 

2-imL  =  0.  (3.31) 

^     ax, 

(3.29)  ensures  that  the  sums  (3.30),  (3.31)  are  finite.  (3.30)  can  then  be  derived  in 

the  following  way.  The  summation  can  be  restricted  to  the  points  in  a  finite  number  m  of 

periods.  2 — ^  'S  the  rate  of  change  of  the  total  area  of  the  Voronoi  polygons  assod- 

ated  with  those  points.  This  total  area  equals  m.  Therefore  (3.30)  follows. 

To  see  (3.31),  let  the  points  2Ci,-X;-  all  move  at  the  same  constant  velocity  u.. 
Then 

E^^  •  u  (3.32) 

is  the  rate  of  change  of  V[2C].  But  V[X.]  clearly  does  not  change  at  all,  P[X.]  is  just 
translated.  Since  u.  is  arbitrary,  (3.31)  follows.  Both  (3.30)  and  (3.31)  can  be  used  to 
represent  the  diagonal  derivatives  in  terms  of  the  off-diagonal  ones.  The  two  representa- 
tions are  different  from  each  other. 

We  turn  to  case  (ii)  now.  We  first  consider  the  derivative  of  V[X.]  in  the  direction 


_lz2L 

I1I-2CI 


.  This  derivative  is 
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jA[X.X] 


as  can  be  seen  from  the  following  figure: 


(3.33) 


area  of  this  sixip: 


Y€A[2C,I]  +  0(e2) 


We  next  consider  the  derivative  in  the  direction  parallel  to  the  face  /*[2C]n^[ll]- 


^2  '^ 


gamed 


Up  to  corrections  of  area  0{i}),  both  the  gained  triangle  and  the  lost  triangle  are  similar 
to  the  triangle  formed  by  2C,Z,Z+ti.  The  scaling  factors  are 


(3.34) 


for  the  gained  triangle  and 


||(2C+I)  -  C[2C,i:]||  +  |a[X,I] 


(3.35) 


for  the  lost  triangle.  Therefore  the  total  area  change  is 


fll2c-rl 


^  (^A[2c,i]-  \\\%^t)-a^xw 


\%-n 


(yA[i,Il+|||(2C+2:)-C[2C,I]||)^ 


(3.36) 
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which  is 

Area  is  lost  moving  in  the  direction  of  the  center  C[2C,2!],  and  area  is  gained  moving 
away  from  Q[^,L]-  From  what  we  have  shown,  it  follows  that 

-^^  =  |A[aC,I]-|j^5^+^^^(^a+I)-£[2C,2:]),  (3.38) 

which  is  equivalent  with  (3.28).  ■ 

We  gave  a  proof  for  proposition  3.7  in  two  space  dimensions.  However,  eq.  (3.28) 
is  correct  in  all  dimensions,  if  C.[X.,X.]  is  defined  as  the  center  of  mass  of  the  face 
between  X.  and  X.-  Peskin's  derivation  of  (3.28)  shows  this;  see  Bttrgers  and  Peskin 
(1985). 

We  notice  that  the  mapping 

(2Ci,---Xv-)   -  (V[2CJ,-,V[20,]) 
is  not  everywhere  twice  continuously  differentiable.   To  see  this,  consider  the  following 

figure. 


•  (2,2) 


With2C:=(l,0)  and  1:=  (1,1  + 8),  we  have  A[2C,2:]=l-8  for  8e[0;l]  here,  endA[2C,I]  =  0 
for  82:1.  Therefore  the  normal  component  of  — 1;™-  is  not  everywhere  differentiable 
(sec  (3.38)). 

To  prove  the  weak  consistency  of  D  with  the  continuous  divergence  operator,  we 
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first  observe  that  the  argument  which  was  used  to  prove  (3.12)  c?-i  be  extended  in  the 
following  way. 

Let  A  be  a  linear  differential  operator.  Assume  that  a  formula  of  the  form 

/A4)(i)4l  =    /(T<J)(i)^  (3.39) 

P  aP 

holds  for  all  smooth  functions  (}>,  where  cr  is  a  linear  differential  operator  on  dP,  possi- 
bly of  order  zero.  In  the  case  where  A  is  the  Laplace  operator,  a<^  =  -^.    Since  cr 

da. 

depends  on  P,  we  also  use  the  notation 

a  =  ap. 

We  discretize  A  by 

V[Ki]L<t>iK,)   :=     Sni*,!],  (3.40) 

where 

I[K^X]   =         /      (T<J)(i)^  +  0(^2).  (3.41) 

If  we  replaced  0{h^)  by  0(h^)  here,  pointwise  consistency  of  first  order  of  L  and  A 
would  follow  immediately.   We  assimie  in  addition  that 

/       <^pgLMi.)'i^  =  -        /       apr^i^(s.)ds  (3.42) 

/[2C,I]  =  -  lUZ]  +  0(h^)  (3.43) 

Proposition  3.8:    (3.39)-(3.43)  together  with  the  non-degeneracy  assumption  (3.11) 
imply  the  weak  consistency  of  L  with  A. 

The  proof  exactly  duplicates  that  of  proposition  3.2.  Analogous  propositions  hold  if 
A  and  L  act  on  vectors  rather  than  scalars,  or  if  the  /»[2C]  are  not  the  Voronoi  polygons 
but  some  different  control  areas  with  mutually  disjoint  interiors,  covering  the  area  of 
interest. 

From  proposition  3.8,  we  now  conclude  the  weak  consistency  of  D. 
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ProposldoD  3.9:    Let  the  non-degeneracy  assumption  (3.11)  be  satisfied.    Then  the 
operator  D  is  weakly  consistent  with  the  continuous  divergence  operator,  i.e. 


J,<^iKi)DiLiXi)V[Ki]  =    /  Hi.)^nis.)di.  +  0(h) 


k  - 


(3.44) 


for  any  smooth  periodic  (scalar  and  vector)  functions  <}>  and  u.. 


Proof: 


I.e. 


n^x]  = 


au=ua  , 


(3.39b) 


y(k(2C)  +  if(I))A[2C,I]^^3^ 


^^f^_ff^  A[;iLX]{\{iL+r)-C.{}LX]) 


(3.40b) 


The  first  summand  in  (3.40b)  originates  from  the  norma]  component  of  (3.38),  the 
second  from  the  tangential  component.  In  the  first  summand,  one  obtains  -u(i)  rather 
than  «(2C)  from  (3.38)  and  (3.31).  To  change  the  sign  is  correct  because  the  integral  over 
BP{lQ  of  the  exterior  unit  normal  is  zero.   We  obtain 


I[KtX]  =        /      u(i)fl^  +  0{h^). 


(3.41b) 


This  relation  would  remain  true  if  we  dropped  the  second  summand  in  (3.40b),  which  is 
itself  of  order  Oihr).  We  thus  obtain  an  alternative  discretization  D;  of  S.--  D^  can  also 
be  derived  by  discretizing  the  divergence  theorem  on  the  Voronoi  polygons  in  the  most 
obvious  way.  See  section  5  for  numerical  resi;lts  concerning  the  comparison  between 
D=D^mdD2. 


Furthermore  we  obtain 


(3.42b) 
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nx.,1.]  =  -  /[IJC]  (exactly).  (3.43b) 


So  (3.44)  follows. 


Proposition  3.10:  D  and  D2  are  pointwise  inconsistent. 

Proof: 

The  inconsistency  of  D  was  first  shown  by  M.  McCracken  (unpublished).  We  use 
the  same  counterexample  as  for  L.  We  consider 

i^^^fc-  (%T)M  VI*   -  i«fe)-.U))  (-S.5)*-  (3-«) 

3«i  24 

The  coefficient  before  - —  in  the  Taylor  expansion  of  (3.45)  is  "-^7*  • 

The  pointwise  inconsistency  of  D2  can  be  shown  with  the  same  counterexample. 

The  coefficient  before  — -  becomes  -  -rrrA^.   ■ 

dx2  425 

We  give  a  second  proof  of  the  inconsistency  of  D .  This  argviment  is  a  reformulation 
of  our  previous  considerations  concerning  the  example  (3.26),  (3.27)  and  shows  a  state- 
ment slightly  more  general  than  proposition  3.10.  As  mentioned  above,  the  area  of  the 
Voronoi  polygon  of  a  marker  at  (0,0)  will,  in  general,  not  be  constant.  This  implies  that 

there  are  fluid  marker  configurations  for  which  Z?u(0,0)9t0,  if  il(s.)=    q 


Take  such  a 


configuration  and  scale  it  by  t.  This  causes  scaling  of  the  gradients  — I'   '   '■  by  t  and 

scaling  of  V[(0,0)]  by  e^.  Because  of  iL(tX.)=inOO,  Du(0,0)  remains  unchanged.  The 
same  argument  shows  the  pointwise  inconsistency  of  any  discrete  divergence  defined  as 
the  relative  rate  of  change  of  the  area  of  a  cell  which  is  scaled  by  e  if  the  marker  confi- 
guration is  scaled  by  e. 
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3.3  Discrete  gradient  operator 

We  define  J2  as  the  negative  adjoint  of  D  with  respect  to  the  discrete  I-  inner  pro- 
duct (3.19).  We  also  define  C:  as  the  negative  adjoint  of  D^  with  respect  to  (3.19). 

Proposition  3.11:  Q,  and  Cj  are  weakly  consistent  with  the  continuous  gradient  operator. 

Proof: 

This  follows  from  the  weak  consistency  of  D,  D2  by  integrating  by  parts.  ■ 

From  (3.30),  it  follows  that  the  kernel  of  Q  contains  the  constants.  (Similarly,  it 
follows  from  (3.31)  that  the  kernel  of  D  contains  the  constant  vector  fields.)  Therefore 
(Z)u,l)  =  0  for  all  u.,  and  L^=Du.  always  has  a  solution  <j>. 

The  kernel  of  G.  may  contain  non-constant  functions.  K,  for  example,  X.y,X\  lie 
on  a  square  grid,  and  if  N  is  the  square  of  an  even  integer,  then  the  kernel  of  Q  is  four- 
dimensional,  for  d  is  then  the  discretization  of  S.  with  central  difference  quotients.  We 
have  been  unable  to  find  a  general  estimate  of  the  dimension  of  the  kernel  of  G.. 
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4.  Solatlon  of  discrete  Helmholtz  eqaadons  on  irregular  (rids 


In  this  section,  we  consider  th?.  fast  numerical  solution  of 


-L«J>  +  c4>  =  /    (c2:0)  . 


(4.1) 


We  use  a  two-level  iteration  which  is  identical  with  the  well-known  multigrid  correction 
cycle,  replacing  the  coarsening  of  the  grid  by  regularization. 

Consider  a  regular  grid  of  the  form 

r*  =  {((i  +  0.5)h,(J  +  0.5)h)  :  ijeZ}  (4.2) 

with  h  =  l/n  and  n=Vv.  i^t  R  be  a  linear  interpolation  operator  which  maps  functions 
defined  on  5v  onto  functions  defined  on  F''.  Let  P  be  a  linear  interpolation  operator  in 
the  opposite  direction. 


•  •       • 

•  •       • 

•  •       • 


We  introduce  the  abbreviations 


A  :=  -L  +  cl  and 


(4.3) 


B  :=  discretization  of  -A  +  cl  on  V^,  using  the  standard  5- point  operator.   (4.4) 


We  consider  the  following  simple  iterative  method  for  (4.1): 


(4.5) 


If  c=0,  S~^  does  not  exist.  B~^r  ii  then  obtained  as  foUows.  Subtract  the  constant  com- 
ponent  from  r,  find  some  solution  v  of  Bv  =  r,  and  subtract  the  constant  component 
from  V.  That  is,  5 "^  is  the  pseudoinverse  of  5. 

If  P  and  R  are  chosen  in  a  simple,  straightforward  way,  the  iteration  (4.5)  con- 
verges very  poorly  if  at  all.    The  reason  is  that  P  and  R  introduce  smoothing,  which 
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prevents  highly  oscillatory  errors  from  converging  fast. 

Therefore  (4.5)  has  to  be  supplemented  by  a  method  which  is  efficient  on  such 
errors.  We  choose  a  suitable  relaxation  method  for  this,  obtaining  the  two-level  cycle 

Step  1:  V]^  relaxation  sweeps 

Step  2:   (4.5)  (4.6) 

Step  3:  Vi  relaxation  sweeps. 

v^  and  V2  are  small  integers,  typically  1  or  2.  To  solve  problems  with  B,  as  required  in 
(4.5),  we  use  the  Fast  Fourier  Transform.  This  will  be  the  only  piece  in  our  fractional 
step  method  which  costs  0(NlogN)  operations  rather  than  0(N)  operations  per  time 
step  Instead  of  the  Fast  Fourier  Transform,  a  very  fast  approximate  solver,  for  example 
a  multigiid  solver,  requiring  only  0(N)  operations  could  be  used.  Such  a  replacement 
often  does  not  have  any  significant  influence  on  the  convergence  factors;  see  Stfiben  and 
Trottenberg  (1982). 

We  have  to  choose  P,  R  and  the  relaxation  scheme.  We  discuss  the  choice  of  P  and 
R  first.  We  list  some  desirable  properties: 

CP:  P  preserves  constants:  P(1)(2C/)  =  1  for  all  j. 

CR:  R  preserves  constants:  R(l)(i)  =  1  for  all  i^r\ 

lOP:  P  preserves  zero  mean  values: 

//24»U)A^  =  0,  then  2:(/'«t»)(2Ci)V[2C]  =  0  . 

lOR:  R  preserves  zero  mean  values: 

//2*(a:*)v^[2C*]  =  0,  then  2:(/?'i>)U)A*  =  0  . 

Here  ^  stands  for        ^       .    lOP  and  lOR  state  that  the  compatibility  conditions  are 
«  xtr'nioap 

preserved  if  c=0.  For  our  two-level  algorithm,  lOR  is  particularly  useful,  since  it  ensures 
solvability  of  the  problems  on  F*  if  c=0. 
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We  also  consider  the  following  related  properties. 
ICP:  P(l)  has  the  discrete  integral  1: 
.     2  (^1)(2^)V[2C]  =  1  . 

ICR:  R(l)  has  the  discrete  integral  1: 
2  (/?l)(i)A^  =  1  . 

IP:  P  preserves  discrete  integrals: 

IR:  R  preserves  discrete  integrals: 

2  (/?4>)(i)A'  =  2  *(2C*)V[2Ct]. 

s  * 

The  following  figtire  shows  trivial  implications. 


CR 

I 


JCR       lOR, 

IR 

We  assume  now  that  P=R',  where  R'  is  the  adjoint  of  R  with  respect  to  the 
discrete  L^-products  on  T''  and  on  S^-,  see  (3.19).  This  means 

(P<^)(Xk)   =   2*^^^^  •t'U)*' and  (4.7) 

{R<t>)U)   =   2c*i'l'(2C*)V[2C*].  (4.8) 

i 

where  the  coefficients  c^^  are  the  same  in  both  equations.  Then  we  obtain  the  additional 
iraplications  shown  in  the  following  figure. 
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CP    ^ 


We  focus  our  attention  on  operators  defined  as  convolutions  of  the  form 

(^4>)(2C)  :=   2  4>(i)8A(iC  -  l.)h\  (4.9) 

(/?4>)U):=    2  *a)8AU  -  i)V[2C].  (4.10) 

where  8^  is  a  spread-out  model  of  the  delta  distribution  with  support  of  linear  size  0{h). 
(4.9)  and  (4.10)  imply  P=R'.  It  is  dear  that  CR  cannot  be  satisfied  by  an  operator  R  of 
this  form,  for  it  is  possible  that  no  point  in  5v  lies  dose  to  x^F'',  in  which  case  R'^is)  is 
0,  no  matter  what  4>  is.  CP,  on  the  other  hand,  is  easily  satisfied,  for  example  with 


^hU)  =  8a(xl)8^(x2),  with 
**■  '  ~   I    0  otherwise. 


(4.11) 
(4.12) 


Notice  that  (4.9),(4.11),(4.12)  is  just  piecewise  bilinear  interpolation.  The  following  fig- 
ure shows  which  of  the  mentioned  properties  the  resulting  operators  R  and  P  have  and 
which  ones  they  do  not  have.  Our  considerations  show  that  no  operator  of  the  form 
(4. 9), (4. 10)  can  have  the  missing  properties. 


Eqs.  (4.9)  •  (4.12)  specify  our  choice  of  P  and  R.  It  remains  to  choose  a  relaxation 
scheme.  The  relaxation  schemes  which  we  shall  describe  are  applied  to  equations  of  the 
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i.e.  we  multiply  the  k-th  equation  by  V[2Ci].  This  has  the  effect  that  the  matrix  is  posi- 
tive definite  and  symmetric  with  respect  to  the  euclidean  inner  product.  It  is  crucial  to 
compensate  for  this  multiplication  in  some  way  when  transferring  residuals  to  the  regular 
grid.  The  same  issue  occurs  when  solving  a  second  order  equation  on  a  regular  grid 
using  a  multigrid  method.  In  that  case,  the  equations  are  often  multiplied  by  the 
squares  of  the  meshwidths  on  all  grids.  A  common  way  of  compensating  for  this  scaling 
is  to  multiply  the  restriction  operator  from  a  grid  with  meshwidth  Aj^  to  a  coarser  grid 

with  meshwidth  A2  ^X  "T-    Similarly,  we  multiply  the  equations  on  the  regular  grid  by 

h^  and  multiply  the  matrix  R  from  the  right  by  the  diagonal  matrix  whose  k-th  diagcnal 

f.2 

entry  is  — ^ — r--  '^^  means  that  we  remove  the  factors  V[^]  before  computing  the  resi- 

duals. 

Using  Gauss-Seidel  relaxation,  we  obtain  a  performance  which  is  disappointing  in 
comparison  with  many  multigrid  methods  on  regular  grids.  The  following  table  shows 
the  reduction  of  the  discrete  Z,^-norms  of  the  residual.  The  order  in  which  the  unknowns 
are  relaxed  is  determined  by  the  numbering  of  the  grid  points,  which  is  random.  In  all 
cases,  the  continuous  problem  being  solved  has  the  solution 

Hi)  =  «n(2Tr(xi  -  2x2))  .  (4.13) 
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Table  4.1.  Residual  reduction  obtained  with  two-level  method  aft^r  1,  3  and  5  itera- 
tions, Helmholtz  equation,  c=10,  V|=2,  V2=l,  Gauss-Seidel  relaxation.  10  random,  uni- 
formly distributed  grids  for  each  N.  Best  results,  worst  results  and  geometric  means 
over  aU  10  experiments. 


N 

worst  case 

best  case 

average 

0.15E-I-00 

0.73E-01 

O.lOE+00 

100 

0.12E-01 

0.17E-C3 

0.12E-02 

0.96E-03 

0.91E-06 

0.22E-04 

0.28E+00 

0.16E+00 

0.22E+00 

400 

0.15E-01 

0.14E-02 

0.67E-02 

0.15E-02 

0.32E-04 

0.38E-03 

0.30E-I-00 

0.21E+00 

0.25E+00 

1600 

0.22E-01 

0.42E-02 

0.98E-02 

0.43E-02 

0.34E-03 

O.llE-02 

The  results  can  be  improved  by  a  modification  of  the  relaxation  method.  This 
modification  is  an  application  of  a  general  rule  due  to  Brandt  (1981)  on  block  relaxation 
in  multigrid  methods.  The  improved  relaxation  scheme  is  defined  as  follows.  Introduce 
the  notation 


When  relaxing  Kt,  determine  the  set  of  all  neighbors  Lot^  for  which 


(4.14) 


(4.15) 


where  8  €  [0;1]  remains  to  be  chosen.  Change  the  values  in  all  these  neighbors  and  in  i* 
simultaneously  so  that  their  equations  become  satisfied.  As  a  result,  some  of  the  ^(Kt) 
may  be  modified  more  than  once  during  one  relaxation  sweep.    We  note  that  our 
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criterion  (4.15)  for  relaxing  points  as  a  block  is  not  equivalent  to  an  analogous  criterion 
based  on  the  distances  between  the  grid  points.  We  have  conducted  experiments  which 
indicate  that  such  a  criterion  leads  to  a  less  efficient  algorithm  than  ours. 

If  8  is  at  least  0.5  and  c^O,  the  blocks  which  are  relaxed  simultaneously  are  at 
most  of  size  2.  To  see  this,  recall  that  the  diagonal  elements  of  L  are  negative,  that  the 
off-diagonal  elements  are  positive  or  zero  with  at  least  3  positive  elements  in  each  row, 
and  that  the  row  sums  of  L  are  zero.  Thus  there  is  at  most  one  neighbor  X.  of  any  point 
Xi  for  which  (4.15)  is  satisfied  when  8^0.5.  With  8=1.0,  we  obtain  the  usual  Gauss- 
Seidel  iteration,  by  a  similar  argument.  We  choose  8=0.5. 

We  repeat  the  above  experiments  with  the  modified  relaxation  scheme. 


Table  4.2.  Experiments  as  in  table  4.1,  with  modified  relaxation  scheme. 


N 

worst  case 

best  case 

average 

0.77E-01 

0.28E-01 

0.48E-01 

100 

0.30E-03 

0.22E-04 

0.72E-04 

0.17E-05 

0.31E-06 

0.65E-06 

0.14E+00 

0.94E-01 

O.llE+00 

400 

0.16E-02 

0.14E-03 

0.40E-03 

0.31E-04 

O.lOE-05 

0.40E-05 

0.16E+00 

0.13E+00 

0.14E+00 

1600 

0.54E-03 

0.27E-03 

0.38E-03 

0.17E-04 

0.58E-05 

0.78E-05 

The  results  for  c=0.0  are  hardly  different.  Here  we  project  the  disoete  right-hand 
side  onto  its  constant-free  part  to  ensure  diat  there  is  a  solution.  It  is  unnecessary  to 
impose  any  condition  such  as  a  prescribed  mean  or  point  value  to  ensure  uniqueness  cf 
the  constant  component  in  the  computed  solution.  The  non-uniqueness  of  the  solution 
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may  simply  be  ignored. 


Table  4.3.  Experiments  as  in  table  4.1,  with  modified  relaxation  scheme,  c=0.0 


N 

worst  case 

best  case 

average 

0.79E-O1 

0.32E-01 

0.53E-01 

100 

0.40E-03 

0.36E-04 

O.lOE-03 

0.26E-05 

0.30E-06 

0.65E-O6 

0.14E+00 

0.98E-01 

0.12E+00 

400 

0.17E-02 

0.17E-03 

0.46E-03 

0.33E-04 

0.13E-05 

0.42E-05 

0.16E+00 

0.13E+00 

0.15E+00 

1600 

0.56E-03 

0.28E-03 

0.39E-03 

0.20E-04 

0.54E-05 

0.78E-05 

We  still  obtain  convergence  factors  which  increase  with  growing  N.  Good  mul- 
tigrid  methods  have  convergence  factors  which  are  bounded  independent  of  the  number 
of  unknowns,  with  a  bound  far  below  1;  see  Stiiben  and  Trottenberg  (1982).  We  do  not 
know  whether  our  method  has  such  a  property.  It  is,  in  any  case,  satisfactory  when  used 
within  our  fractional  step  method,  see  section  6. 
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5.  Projection  of  a  vector  field  onto  the  kernel  of  the  discrete  divergence  operator 

The  last  tool  needed  for  our  fractional  step  method  is  an  algorithm  which  projects  a 
given  vector  field  u.  on  5v  orthogonally  onto  the  kernel  of  the  discrete  divergence  opera- 
tor D.  This  means  that  we  want  to  find  Pu  and  q  such  that 

U  =  Vu  +  aq   and  (5.1) 

D?u.  =  0  .  (5.2) 

For  this  purpose,  we  have  to  solve 

DQq  =  Dii  .  (5.3) 

DQ  is  singular,  but  Du.  satisfies  the  compatibility  condition,  i.e.  Du  is  orthogonal  to  the 
kernel  of  DQ.  is  the  inner  product  (3.19),  and  Qjq  is  unique. 

On  a  square  grid  with  meshwidth  h,  DQ  is  the  standard  5-point  discretization  of  the 
Laplace  operator  with  meshwidth  2h.  Simple  relaxation  schemes  have  no  smoothing 
effect  for  this  discretization,  as  can  immediately  be  verified  by  Fourier  analysis.  We 
ignore  the  fact  that  the  system  can  be  decoupled  into  4  smaller  systems  which  can  be 
solved  easily,  since  this  will  not  be  the  case  on  irregular  grids.  The  difficulty  is  related  to 
the  fact  that  there  are  highly  oscillatory  grid  functions  in  the  kernel  of  the  finite  differ- 
ence operator,  and  such  components  are  not  affected  at  all  by  relaxation  methods.  By  a 
continuity  argument,  the  error  in  a  nearby  Fourier  mode  is  damped  slowly.  Problems  of 
this  kind  have  been  investigated  much  more  generally  by  Brandt  and  Dinar  (1979). 

One  might  consider  solving  (5.3)  by  a  conjugate  gradient  iteration  with  precondi- 
tioning. On  a  square  mesh,  it  is  useless  to  precondition  DQ  with  L.  This  can  immedi- 
ately be  shown  by  discrete  Fourier  analysis.  Even  if  a  good  preconditioner  for  DQ.  could 
be  found,  the  fact  that  we  do  not  know  much  about  the  kernel  of  DQ  ,  i.e.  the  kernel  of 
G.J  would  present  a  problem.  A  non-trivial  null-space  does  not  affect  the  performance 
of  the  conjugate  gradient  method  in  exact  arithmetic.  In  floating  point  arithmetic,  how- 
ever, the  method  often  does  not  converge  when  applied  to  a  positive  semidefinite  matrix 
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with  a  non-trivial  null-space  unless  the  computed  residuals  are  projected  into  the  orthog- 
onal complement  of  the  null-space. 

We  therefore  replace  the  discrete  projection  operator 

/  -  a{Da)-'D  (5.4) 

by 

/  -  Gl^-^D  .  (5.5) 

DQ  and  L  are  singular.  Nevertheless  Q.{DQ.)~^D  and  QU'^D  have  obvious  well-defined 
meanings. 

(5.5)  is  not  a  projection  operator.  In  fact,  the  discrete  L*-norm 

11/  -  GZ.-^D||  =  p(/  -  QL-'D)  (5.6) 

is  often,  i.e.  on  many  fluid  marker  configurations,  larger  than  1.  The  eigenvalues  with 
largest  absolute  value  must  clearly  be  smaller  than  —  1  if  that  is  the  case. 

Proposition  5.1:  If  all  eigenvalues  of  I—Ql.~^D  lie  in  (-1;1],  then 


lim  (/  -  Gl'^oy  =  I  -  a(Da)~^D  .  (5.7) 

J-" 


Proof: 


£  :=  lim  (/  -  QJ.~^Dy(u.)  differs  from  u  only  by  a  discrete  gradient.  Therefore  it 

suffices  to  prove  that  Dii=0.  Qearly,  (I-QJ^'^D)^  =  i,  therefore  Ql~^D^=0.  Using 
Q.=  -D',  we  conclude  that  L~^D^=0,  therefore  D^=0.  As  before,  L~^  has  an  obvious 
meaning.   ■ 

We  give  a  different  interpretation  of  iI-QJ^~^D)lu..  To  solve  (5.3),  solve  Lq^=Du. 
first  and  perform  7- 1  defect  correction  iteration  steps  (sec  Stetter,  1978),  using  L  on  the 
left-hand  side  and  DQ,  on  the  right-hand  side,  obtaining  an  approximation  qj.  Then 

ii-aq,  =  (I-Ql-'Dyii . 
Since  (5.5)  is  a  discretization  of  the  continuous  projection  operator,  it  is  to  be  expected 

that  at  least  those  eigenvalues  which  correspond  to  smooth  eigenfunctioss  are  larger  than 
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-1.  This  motivates  the  following  modification  of  (5.5): 

/  -  C(/  -h  uiiyL-^D  ,  (5.8) 

with  <D~l/p(I,)  but  a)Sl/p(I,),  k  integer.   In  our  experiments,  we  always  take 

w  :=  T  .  (5.9) 

P 

where  p  is  the  maximum  row  sum  norm  of  the  matrix  L,  i.e. 


u>  = 


max  2 


'i-i 


*  niUii»K,\\i^-^\ 


(5.10) 


We  show  that  o)  is  at  most  of  order  0{h^).  Using  the  estimates 

u-n^2h  (5.11) 


2A[i,I]a:2V:|^V^  (5.12) 

1*\ 


in  (5.10),  we  obtain 


AminVVgl 

"^    2v;-    •  (^-^^^ 

/-l-a)L  is  therefore  pointwise  consistent  with  the  identity,  since  L^  is  bounded  for  h-^  if 
({)  is  a  fixed  smooth  periodic  function;  see  proposition  3.4. 

For  fixed  k,  (5.8)  can  still  be  considered  a  discretization  of  the  continuous  projec- 
tion operator. 

Proposition  5.2:    For  sufficiently  large  k,  the  powers  of  (5.8)  converge  to  the  exact 
discrete  projection  operator  I  -  Q.{DQ,)~^D. 

Proof: 

The  proof  that  the  limit,  if  it  exists,  is  /  -  0.(0(2)' ^D  is  a  duplicate  of  the  proof  of 
(5.7).  To  see  that  the  limit  of  the  powers  of  (5.8)  exists,  notice  that 
((/-C(/-t-toL)*L"^D)  is  symmetric  in  the  discrete  L^-product  (3.19).  We  show  that  the 
Rayleigh  quotient 
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((i-a(i+i^Ln-'D)uM) 

lies  in  (-1;1]  for  large  k.  For  this  purpose,  we  rewrite  (5.14)  as 


(        ^  (5.14) 


1    -H     ((/^coL)*L-^D|t^M)  (5.15) 

Since  (/-i-wZ,)*!"^  is  symmetric  and  negative  semidefmite  with  respect  to  (3.19), 
(5.15)  is  si.  For  large  k,  (5.15)  is  >-l,  because  those  eigenvalues  of  (Z+uZ,)*  which 
are  different  from  the  simple  eigenvalue  1  tend  to  0  for  k-^.  ■ 

This  leads  to  the  following  algorithm,  which  chooses  the  smallest  possible  k  itself: 

Algorithm  5.1: 

Given  u.,  generate  uP,  n^,  u^,---  with  i*/  -  Pu  as  follows: 

14°  :=  u  . 

Given  it!  ,  define 

*  :=  0;    ii/'°  :=  u/,  <f'^  :=  L-^Du/'^;    u!'^  :=  tt/>°  -  aq>'K 

H'/i//*  ||lt/-**MI>  lltt'il: 

In  our  experiments,  this  algorithm  usually  requires  k=0  when  computing  u^  from  n^, 
k>0  for  j>l,  sometimes  unfortunately  k  »  0. 

Notice  that  |iit''^^||s|ltt'||.  This  property  ensures  the  linear  stability  of  the  fractional 
step  method  introduced  in  section  6,  even  if  only  a  few  iterations  of  the  projection  algo- 
rithm 5.1  are  used. 

One  can  interpret  the  factor  (I+inLY  in  the  following  way.  Application  of  /-fuL 
to  a  fimction  q^  is  equivalent  to  performing  one  modified  Jacobi  relaxation  sweep  for  the 
system  Lq  =  0,  with  the  initial  guess  ^q.  We  expect  such  a  sweep  to  damp  highly  oscilla- 
tory components  of  q^  and  leave  smooth  components  essentially  unchanged.  It  might 
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therefore  be  useful  to  replace  the  factor  I+oiL  by  one  sweep  of  a  more  efficient  smooth- 
ing method.  Numerical  results  indicate,  however,  that  the  relaxation  method  used  in  sec- 
tion 4  is  significantly  inferior  to  modified  Jacobi  relaxation,  i.e.  to  the  smoothing  factor 
/-t-uL,  in  the  present  context. 

We  give  a  nimierical  example.  Consider 

u(s)  =  (jin(2iTXi)cos(2iTX2)  ,  0)  .  (5.16) 

The  divergence-free  part  of  u,  is 

0.5(jtn(2iTaC|)cos(2'iTX2)  ,  -cos(2iTxJsin(2iTj:2))  .  (5.17) 

We  compute  the  divergence- free  part  of  u.  numericaUy,  compute  the  discrete  L^-norms  of 
the  discretization  errors  in  the  two  components,  and  compute  the  square  root  of  the  sum 
of  the  squares  of  these  two  numbers.  This  results  in  the  following  table. 
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Table  5.1.    Discretization  errors  obtained  witb  algorithm  5.1  for  example  (5.16).    20 
random,  uniformly  distributed  grids  for  eadi  N. 


N 

-iteration 

worst  case 

best  case 

N 

iteration 

worst  case 

best  case 

1 

0.30 

0.22 

1 

0.081 

0.070 

2 

0.30 

0.20 

2 

0.078 

0.061 

3 

0.28 

0.20 

3 

0.083 

0.063 

25 

4 

0.29 

0.19 

400 

4 

0.081 

0.064 

5 

0.29 

0.20 

5 

0.086 

0.064 

10 

0.29 

0.20 

10 

0.084 

0.062 

15 

0.30 

0.20 

15 

0.089 

0.062 

1 

0.17 

0.13 

1 

0.041 

0.035 

2 

0.16 

0.11 

2 

0.040 

0.031 

3 

0.16 

0.11 

3 

0.033 

100 

4 

0.16 

0.11 

1600 

4 

0.034 

5 

0.16 

0.11 

5 

0.035 

10 

0.16 

0.11 

10 

0.037 

15 

0.16 

0.11 

15 

0.037 

"•"  means  that  k=50  is  not  sufficient.  In  none  of  the  experiments  summarized  in  the 
table,  k  is  larger  than  0  in  the  first  iteration  or  larger  than  2  in  the  second  iteration. 

One  iteration  generates  an  approximation  to  the  continuous  solution  which  is  as 
good  as  or  even  better  than  the  approximation  obtained  after  many  iterations.  What  we 
have  gained  by  replacing  (5.5)  with  (5.8)  are  guaranteed  stability  and  convergence  to  the 
solution  of  the  discrete  projection  {problem,  not  higher  accuracy. 

It  is  surprising  that  the  truncation  error  does  decrease  with  growing  N  here,  even 
though  D  is  not  pointwise  consistent  with  the  divergence  operator.  We  also  observe  that 
the  variance  of  the  trimcation  error  seems  to  decrease  with  growing  N.    These  two 
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statements  are  based  on  numerical  results.  We  have  not  been  able  to  prove  them. 

The  operator  D^  introduced  in  section  3,  i.e.  the  operator  which  is  obtained  by 
discretizing  the  divergence  theorem  on  the  Voronoi  polygons  in  the  most  obvious  way, 
appears  to  give  far  worse  results  here.  We  illustrate  this  by  recomputing  a  part  of  the 
previous  table,  now  using  Di  rather  than  D=D^  and  jCZ2»  the  negative  adjoint  of  D2, 
rather  than  C=l2i- 

Table  5.2.  Experiments  as  in  table  4.1,  with  D, (2  replaced  by  02,0.2- 


N 

iteration 

worst  case 

best  case 

N 

iteration 

worst  case 

best  case 

1 

0.29 

0.22 

1 

0.11 

0.095 

25 

2 

0.27 

0.18 

400 

2 

0.098 

0.085 

1 

0.16 

0.12 

1 

0.094 

0.086 

100 

2 

0.14 

0.11 

1600 

2 

0.083 

0.076 
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6.  A  b-actional  step  method  for  the  Navier-Stokes  eqaatioiu 

We  consider  the  foUowing  fractional  step  method  for  the  Navier-Stokes  equations. 

^"  +  1  =  Pi"*!  (6.1) 

Here  L"  denotes  the  discrete  Lapladan  on  i'v,  the  set  5v  at  time  nit.  The  points  in  that 
set  are  the  components  of  the  vector  2C" .  P"  is  the  orthogonal  projection  onto  the  kernel 
of  D",  where  D"  is  the  discrete  divergence  operator  on  S^,-.  The  pressure  is  treated  in 
(6.1)  as  in  the  fractional  step  methods  introduced  and  analyzed  by  Chorin  (1968,1969) 
and  Temam  (1969).  The  nonlinear  term  is  absent  in  this  Lagrangian  model. 

Writing  PV*^  =:  i"*^-  C<7''*^  we  obtain 

'^   ^^'^    -  vL-i"^^  +  ^^^      =  r*'    and  (6.2) 

^nj^/i  +  l  =  0  . 

Therefore  p""*"^  :=  q"'*'^/lt  is  an  approximation  of  the  pressxire.  The  formulas  (6.2)  with 
L''u.'''*'^  instead  of  L"u"'*'^,  together  with  the  third  equation  of  (6.1),  define  the  method 
introduced  by  Peskin  (1979). 

(6.1)  is  unconditionally  linearly  stable  in  the  sense  that 

llu''*^||(„)^Ilu''||(„)  +  Af|lr*^||(„)  (6.3) 

where  ||"||(„)  denotes  the  discrete  L^-norm  on  5Jv-.  This  follows  from  proposition  3.5. 
Using  the  projection  algorithm  introduced  in  section  5,  stability  is  guaranteed  even 
though  we  do  no  apply  the  operator  P"  exactly. 

We  also  consider  the  following  modification  of  (6.1). 
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Notice  that  this  method  requires  two,  not  three  applications  of  P"  per  time  step. 
Writing,  as  before,  P"u"*^  =:  il"'*'^-QPq''*^,  and  writing 

r"(u"  +  ^t£'*^)  =:u."  +  \t£'-^^-a"r"-*-\ 
we  obtain 

JL_J^  -  vL'-ii"^^  +  a"i^—jf—-)  =ir^^    and  (6.2b) 


^«^«+i  _  0 


Therefore/?"*^  :=  -» — is  an  approximation  of  the  pressure. 

Those  of  the  following  nimierical  results  for  which  the  method  is  not  specified  are 
obtained  with  method  (6.1). 

Teat  problem  1: 

We  construct  an  exterior  force  density  such  that  the  flow  becomes 

«i(i,0  =    Jin(Yr)coj(2iTJcJain(2iTJ2) 

"lU'O  =  -sin(^t)sm(2'ttx^cos(2vx2)  (6.4) 

p(i.O  ~  sin(^t)cos(2'ttx^cos(2'nx2). 
We  measure  the  discrete  Z,^-norms  of  the  errors  in  velocity  and  pressure  at  time  1. 
In  every  time  step,  we  use  only  one  step  of  algorithm  5.1.  Initially,  the  ^  are  at  the 
positions  ((k^  +  0.5)h,(k2  +  O.S)h),  Ostj^sV^-1,  O^skj^^-l.  We  expect  that  the 
solution  of  the  scalar  Helmholtz  and  Poisson  problems  up  to  rounding  error  accuracy  is 
necessary  only  at  time  0.  At  a  later  time  step,  the  values  from  the  previous  time  step 
should  provide  an  excellent  initial  guess,  and  little  work  should  be  required  to  improve 
this  guess  such  that  the  truncation  error  level  is  reached.  We  use  only  one  two-level 
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cyde  per  scalar  problem,  with  v^=2,  V2=l.  The  following  tables  show  results  of  oxir 
experiments.  We  use  the  abbreviations  /ly  :=      r—,  e^  :=  error  in  u^,  e-^  :=  error  in  u^ 

and  e^  -:=_  error  in  p.  Since  the  average  of  the  exact  pressure  is  zero,  we  project  the 
computed  approximation  for  the  pressure  onto  its  constant  free  part  with  respect  to  the 
discrete  L^-product  before  measuring  e^. 

Table  6.1.  Errors  in  velocity  and  pressure,  test  problem  1,  ^t=hf{,  v=1.0,  starting  on 
square  grids. 


N 

*i 

«2 

'p 

100 

0.35 

0.35 

0.11 

400 

0.12 

0.12 

0.38 

1600 

0.063 

0.064 

0.33 

From  these  results,  it  seems  unlikely  that  the  approximation  for  the  pressure  con- 
verges to  the  exact  pressure  at  all  in  this  case. 

We  repeat  the  experiments  with  lower  viscosity. 

Table  6.2.   Errors  in  velocity  and  pressure,  test  problem  1,  At=h^;  v=0.1,  starting  on 
square  grids. 


N 

«i 

«2 

^P 

100 

400 

1600 

0.35 
0.12 
0.064 

0.35 
0.12 
0.065 

0.11 
0.12 
0.082 

The  following  table  contains  results  obtained  with  random  initial  grids.   As  before, 
we  use  one  two-level  cycle  per  scalar  problem.  Notice  that  the  errors  do  not  differ  very 
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much  from  those  displayed  in  table  6.2. 

Table  6.3.   Errors  in  velocity  and  pressure,  test  problem  1,  Ar=Av.  v=0.1,  starting  on 
20  imiformly  distributed,  random  grids. 


N 

<i 

*2 

'p 

worst 

best 

average 

worst 

best 

average 

worst 

best 

average 

100 

0.51 

0.31 

0.38 

0.48 

0.29 

0.36 

0.37 

0.20 

0.27 

400 

0.17 

0.13 

0.15 

0.18 

0.13 

0.15 

0.14 

0.086 

0.12 

1600 

0.064 

0.059 

0.061 

0.062 

0.059 

0.061 

0.076 

0.069 

0.072 

As  in  section  5,  we  observe  that  the  variance  of  the  error  appears  to  decrease  for 
increasing  N. 

We  repeat  some  of  the  above  experiments,  now  solving  all  scalar  problems  up  to 
truncation  error  accuracy,  but  still  applying  only  one  step  of  algorithm  5.1.  The  results 
confirm  that  one  two-level  cycle  per  scalar  problem  is  sufficient: 


Table  6.4.  Experiments  as  in  table  6.2,  solving  all  scalar  problems  up  to  roimding  error 
accuracy. 


N 

«i 

«2 

'? 

100 

0.35 

0.35 

0.11 

400 

0.12 

0.12 

0.12 

1600 

0.067 

0.068 

0.086 

The  slight  difference  between  the  errors  in  the  two  velocity  components  for  N=  1600  is 
caused  by  rounding.  In  exact  arithmetic,  tfi=ff2- 
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Using  smaller  viscosity  coefficients,  for  example  v=10"*,  we  still  obtain  reasonable 
results  for  this  problem.  However,  some  of  oixr  results  for  test  pjroblem  4  below  indicate 
that  the  method  cannot  generally  be  used  for  such  viscosity  coefficients. 

Test  problem  2: 

We  consider  a  shear  flow; 

"2(i.O  =    0  (6.6) 

PU,t)  =    0 

PropositloD  6.1:  If  the  initial  marker  configuration  is  a  rectangular  grid  with 
meshwidths  Ax^  and  Ax^, 

AxjayAxi,  (6.7) 

and  if  the  initial  velocity  and  the  exterior  force  are  independent  of  x^  and  parallel  to  the 
x^-axis,  then  the  discrete  approximation  for  the  second  velocity  component  obtained  by 
(6.1)  or  (6.1b)  is  0. 

Proof: 

The  statement  would  be  obvious  if  P"  were  replaced  by  the  identity  operator. 
Therefore,  it  suffices  to  prove  the  following  statement. 

Assume  that 

{2Ci.-J^}=  {((«  +  y)Ax,+8^,0  +  y)Axj)  :.  =  0,--,-^-l,  ;  =  0.  •■•,-^-l}(6.8) 
with  numbers 

8;€[-i-Axi,lAxJ  (6.9) 


independent  of  i,  and  assume  that  (6.7)  holds.   Consider  a  vector  field 


(6.10) 
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Then 


DyiiKi)  =  0 


(6.^  } 


for  all  k.  - 


To  prove  thit,  consider  two  rows  of  fluid  markers  and  the  associated  Voronoi 
diagram  in  an  infinite  (periodic)  horizontal  strip,  as  indicated  in  the  following  figure. 


«      • 

•      • 

• 

• 

•    • 
I 

• 

•     • 

. 

• 

[  r 

•  •  •  i « 1  • 

! 

All  Voronoi  polygons  have  the  same  shape  and  therefore  the  same  area. 

Now  we  consider  the  case  of  more  than  two  rows  of  markers.  (6.7)  ensures  that 
the  Voronoi  polygons  in  row  j+1  are  not  adjacent  to  the  Voronoi  polygons  in  rowj  — 1. 
The  Voronoi  polygons  in  row  j  are  similar  to  each  other,  and  the  total  area  which  they 
occupy  within  one  period  is  constant.  Therefore  each  individual  polygon  has  a  constant 
area,  which  implies  (6.11).  ■ 

We  start  the  calculation  on  random  grids.  We  obtain  the  following  results. 
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Table  6.5.   Erron  ia  relodty  and  pressure,  test  problem  2,  Ar=A,^-,  v=1.0,  startiijg  on 
10  random,  mnfotmly  distnbuted  grids. 


N 

"    -                «i 

*t                       ! 

*? 

worst 

best 

average 

won-t 

best 

average 

worst 

best 

average 

100 

400 

1600 

0.075 
0.035 
0.016 

0.045 
0.025 
0.014 

0.060 
0.030 
0.015 

0.071 
0.035 
0.015 

0.041 
0.024 
0.013 

0.056 
0.029 
0.014 

0.061 
0.033 
0.015 

0.041 
0.022 
0.013 

0.053 
0.027 
0.024 

Test  problem  3: 

The  solutions  u  of  problems  1  and  2  have  the  spedal  property  that  wYjl  is  a  gra- 
dient. Therefore  we  test  the  method  on  another  problem  with  a  simple  prescribed  solu- 
tion which  does  not  have  this  property: 

Mi(i,r)  =  sin2(2irj2) 

«2U.O  =  cos(2ttxJ  (6.12) 

p{i.,i)  =  sin^(2'iTxJ  +  sin2(2iTac2) 

At  time  0,  we  solve  the  scalar  problems  up  to  rounding  error.   At  later  time  steps,  we 
solve  them  approximately,  using  one  two-level  cycle.  We  obtain  the  following  results. 


Table  6.6.   Errors  in  velocity  and  pressure,  test  problem  3,  Ar=Ay,  v=0.1,  starting  on 
square  grids. 


N 

«i 

*2 

'P 

100 

0.31 

0.25 

0.47 

400 

0.14 

0.12 

0.27 

1600 

0.070 

0.058 

0.16 
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Test  problem  4: 

We  consider  an  initial  vortex  blob  centered  at  (0.5,0.5)  and  compute  the  flow 
which  results  in  the  absence  of  any  exterior  force.  To  obtain  the  initial  velocity  Held,  we 
solve  the  discrete  Poisson  problem  with  periodic  boundary  conditions  and  the  right-hand 
side 

'    1    (l+cos(^I^^i^))(l+cos(^%;^))  +  c  if  lx,-0.5|50.2  and  lx2-0.5|^0.2 


0.16'  '       0.2       "^  ^       0.2 

,       .  (6.13), 

c  otherwise 


where  the  constant  c  is  chosen  such  that  the  discrete  compatibility  condition  is  satisfied. 
From  the  resulting  discrete  stream  function  iji'',  we  obtain  the  velocity  using  central  dif- 
ferencing. 

Our  numerical  results  for  this  example  show  that  the  centrifugal  force  tends  to  push 
the  markers  away  from  the  center,  leaving  a  region  in  which  there  are  no  markers  at  all. 
To  counteract  this  tendency,  we  use  the  modified  method  (6.1b),  in  which  the  velocity 
field  is  projected  onto  its  divergence  free  component  not  only  before  but  also  after  mov- 
ing the  markers.  Again  we  use  only  one  step  of  algorithm  5.1  in  each  projection. 

The  following  figures  show  computed  flow  fields  for  N=400,  v=0.1,  Ar=0.025,  at 
time  t  =  0.0,  0.1,  0.2,  0.3.  Ar  was  0.025  here.  Due  to  viscosity,  the  velocity  decreases 
rapidly.  The  maximum  computed  speed  is  0.89  at  time  0.0  and  0.14  at  time  0.3.  This 
decay  is  not  visible  in  our  plots  since  we  normalize  the  lengths  of  the  arrows  such  that 
the  maximum  length  is  the  same  for  each  plot.  The  plots  obtained  at  later  time  do  not 
differ  much  from  the  one  at  time  0.3. 
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We  repeat  this  computation  with  v=0.01  and  obtain  the  following  fields  at  the 
same  times: 
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We  perform  the  following  convergence  test  We  perform  the  computations  with 
N=  100,400  and  1600,  each  time  starting  on 

Notice  that  the  initial  marker  sets  for  N=400,  N=1600  contain  the  initial  marker  set  for 
N=100  in  this  case.    Let  2Ci'*,"JCi§S  be  the  marker  positions  for  N=100  at  time  0.3. 


•71- 

Let  2C/*^  and   Xi^^  be  the  positions  of  the  same  markers  at  time  0.3  computed  with 
N=400  and  N=1600.  We  measure 


100 

/-I 


•100,400 


and 


100 


100 

/-I 


For  a  first  order  convergent  method,  one  would  expect 


^400,1600  ' 2 — • 

The  computed  values  are  even  slightly  more  favorable: 


Table  6.7.  Test  problem  4,  «ioo.400  and  -r^ocieoo.  ^'=^ 


V 

'100.400 

'400.1600 

0.1 

0.0090 

0.0040 

0.01 

0.065 

0.018 

0.001 

0.12 

0.043 

0.0001 

0.17 

0.086 

Thus  convergence  for  N-<^  is  apparently  obtained  even  for  very  small  viscosity 
coefficients.  Notice,  however,  that  the  size  of  the  error  increases  substantially  for  v-^ 
when  N  is  fixed. 
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7.  Elastic  booiMUrles  Immened  In  the  flnld 

In  this  section,  we  report  on  some  preliminary  numerical  experiments  with 
immersed  boundaries.  The  mathematical  model  and  the  numerical  method  which  we  use 
here  are  due  to  Peskin  (1977);  see  also  Pesldn  and  McQueen  (1980). 

We  first  present  the  mathematical  model.  Consider  ilow  in  a  region  Q  with  an  elas- 
tic, massless  boundary  immersed  in  the  fluid.  Let  2C(«,0  ^  ^  representation  of  this 
boundary  at  time  t  such  that,  for  fixed  s,  i{s,t)  is  the  trajectory  of  a  fixed  boundary 
point. 

The  boundary  exerts  a  force  on  the  fluid  which  can  be  determined  in  the  following 
way.  Consider  an  arc  s€(a;b)  of  the  boundary.  The  mass  of  this  arc  is  zero,  therefore 
its  momentimi  and  the  rate  of  change  of  its  momentum  are  zero.  The  sum  of  the  forces 
exerted  on  it  must  therefore  be  zero.  The  fluid  exerts  a  force  on  the  arc,  which  we 
assimie  to  be  of  the  form 

-Uii,t)d,.  (7.1) 

Therefore,  the  force  exerted  by  the  arc  on  the  fluid  is 

Ui',t)ds  (7.2) 

A  second  force  is  exerted  on  the  arc  by  the  boimdary  itself.  It  is  eqxial  to 

Tit  (7.3) 

where  T  is  the  tension,  defined  by 

r  :=  5(11-^^^11  -  1)  ,    and  (7.4) 

1  :=      J',    .  (7.5) 

S  is  a  material  constant,  a  stiffness  coefflcient.  We  conclude 

b 

Jl(s,t)d3  =  Tzt  (7.6) 
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Fonnally,  our  mathcmatica]  problem  is  then 

u,  +  (mDm  -  vAi4  +  Sp   =  £U,r)  (7.7) 

Sit  =   0 

uU.O)   =  Q 

where  E(s.,t)  is,  for  fixed  c,  the  disti4bution  whose  action  on  vector  valued  test  functions 
is  defined  by 

iE,i^)   :=   J^^MU',t))d3.  (7.8) 

We  have  assumed  here  that  the  density  p  eqiials  1.  A  change  in  p  has  the  effect  of 

c 

changing  5  in  this  model.  Only  —  is  relevant. 

Next  we  describe  the  numerical  method.  We  consider  a  case  in  which  the  boundary 
is  a  dosed  simple  curve.  It  is  represented  by  M  fluid  markers  2Ci."*  JCv^-  The  remaining 
N-M  markers  represent  the  fluid.  We  number  the  particles  such  that  the  j-th  marker  is 
next  to  the  (j-l)-st  and  (j+l)-st  markers  in  the  circular  diain  of  markers  representing 
the  boundary.  In  this  section,  index  calculations  are  carried  out  modulo  M. 

We  number  the  boundary  markers  first  for  notational  convenience  only.  In  our 
code,  the  boundary  markers  are  actually  numbered  last.  This  has  slight  advantages  for 
the  construction  of  Voronoi  diagrams,  since  the  Voronoi  diagrams  associated  with  the 
boundary  markers  only  are  degenerate. 

For  j^{l;--]M},  the  force  on  P[X.j]  is  defined  by 

^l      A.  'J  |lx,.,-2(,||  ^  1      A.  'jlH^-r^r      ^^-^^ 

which  discretizes  Til^.  Aj  is  an  input  parameter,  the  resting  length  of  the  springs  which 
form  the  boundary.  Aj  determines  the  possible  resting  positions  of  the  boundary.  A 
corresponding  input  to  the  continuous  model  is  the  parametrization  of  the  initial  curve 
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2Co(j). 

Following  Peskm  (1977),  we  compute  the  forces  exerted  by  the  boundary  on  the 
fluid  inthe  following  implicit  way.  Given  marker  positions  2C/  and  velocities  Uj  at  time 
nAf  (j=l,...,N),  we  first  compute  approximate  new  positions  2Ci,-"  JC,v  of  the  boundary 
markers.  These  approximate  new  positions  are  only  used  to  compute  the  force  exerted 
by  the  boundary  on  the  fluid.  The  actual  new  positions  are  then  obtained  by  advancing 
X.j^llj  by  one  time  step  in  exactly  the  same  way  as  in  section  6,  using  the  computed  force 
exerted  by  the  boundary  as  the  right-hand  side. 

We  compute  the  auxiliary  positions  i* ,  •  •  •  JQ  ^^om 


2C/*  =  2Cy  +  Arty  + 


JML 


u;.i-m 


Aj 


u 


-1 


Il2g%x-2c;i 


+  s 


As 


-] 


^-^    ^ 


iix;-i-i;i 


(7.10) 


As  in  the  continuous  case,  the  density  is  assumed  to  be  1  here,  and  a  diange  in  density 
can  be  interpreted  as  a  change  in  the  stiffness  S. 

We  introduce  the  following  simplifying  notation: 


^''     [v[^"]^s\    ' 

(7.11) 

2C  :=  (2Ci,""»2C,vf)» 

(7.12) 

y(Vi-2^)^  -  ^s\\Kj^,-lj\\ 

(7.13) 

X^  :=  iQ+AtU^. 

(7.14) 

J 


(7.10)  can  then  be  rewritten  as 

2C;  =  2^-5y2Si?(r).  (7.15) 

As  described  so  far,  the  disaete  model  describes  the  boundary  as  a  circular  chain 

of  springs  which  have  a  resting  length  Aj  and  resist  compression  as  well  as  extension. 

Peskin  (1977)  showed  that  the  solution  of  (7.15)  becomes  much  simpler  if  the  resistance 

against  compression  is  dropped.   This  means  that  we  modify  the  definition  of  E  in  the 
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following  way: 


EU)  :=  2 

J 


(7.16) 
0  otherwise. 


E  is  then  once  but  not  twice  continuously  differentiable. 

With  this  modification,  the  solvability  of  (7.15)  can  be  established  and  a  solution 
can  be  found  in  the  following  way.  Let 

i:°:=(«r^2C?.-,iw^2C^)  (7.17) 

and  consider 

EU)  :=  yIll-I<'|P  +  £(5ii::.-,UXv).  (7-18) 

a  function  of  Y.=  iZi,--,Zu)-   E  has  at  least  one  local  minimum,  since  E^O.  Let  Ji  be  a 
local  minimum  of  E.  Then 

1Q  =  ^-  €y(I^)(5:i:i.-,ivXv)-  (7.19) 

Let 

Xj  :=  ^jlj,  (7.20) 

then 

2(;/  =  i?-5;Y^(2Ci,   -JCv).  (7.21) 

An  analogous  change  of  coordinates  can  be  used,  more  generally,  for  systems  of  the 
form 

X.  =  ie-^^(^)  (7-22) 

with  a  positive  deflnite  matrix  A  and  a  function  E  boimded  from  below. 
We  introduce  the  notation 

hj  :=  ^y^iI/^i-Syr,.  (7.23) 

Then  E(X.)  can  be  written  as 
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0  otherwise. 


(7.24) 


We  use  Newton's  method  with  line  search  to  minimi/ft  (7.24).  For  this,  we  com- 
pute the  gradient  and  the  Hessian  matrix  of  (7.24).  The  gradient  is 


The  Hessian  matrix  is 


-max( 


(llS^||-A.,0)-i^+max(|lfi^_J|-A^,0)-i||:^ 


I+(^,^jC,j),i=l,...^fJ=l,...J^, 
where  the  Cjj  are  2x2  blocks  given  by  the  following  formulae. 


(7.25) 


(7.26) 


C,y  =  0  if   |(«-;)(modM)|a2 


^ij  ~     ^ij-i      Q;+i 


Q+ij  ~  Qj+ 


C/ ,+1  =  -/+ 


Aj 


llfi^lP 


(7.27) 
(7.28) 
(7.29) 

(7.30) 


with  the  term  [•••]  in  (7.30)  replaced  by  0  if  |lfij|I<Aj.  Notice  that  all  C/j  are  symmetric. 

Pesldn  (1977)  proved  that  the  Hessian  is  positive  definite.   We  present  this  proof  in 
our  case. 

Proposition  7.1:  The  matrix  (7.26)  is  positive  defmite. 

Proof: 

First,  we  rewrite  the  Hessian  as 


/  +   diag(5,)  ■  C  •  diag(?,), 


(7.31) 


where 


C   :=    (C, ,) 


(7.32) 


-77 


It  suffices  to  prove  that  C  is  positive  semidefinite.  We  first  show  that  the  matrices  C,  ,+i 
are  negative  semidefinite.  Cij+i  differs  from  -/  by  a  singxilar  summand,  so  one  of  its 
eigenvalues  is  —  1.  The  sum  of  the  eigenvalues  of  Cf^^^  is 


trace(C,  ,^J  =  -2  + 


^s 


-2  otherwise. 


(7.33) 


Therefore  the  second  eigenvalue  must  be  nonpositivc. 

We  conclude  from  this  that  C  is  positive  semidefinite.  Let  di,-,d,i4  be  vectors  of 
length  2.  We  have  to  show  that 


/  /  / 


(7.34) 

is  nonnegative.  Subscript  calculations  are  always  carried  out  modulo  M.   Using  (7.28), 
we  rewrite  (7.34)  as 

-  '2d[C,j^,di  -  ^d[C,j.,-d,  +  Sir-C,;^:-^,^  +  SiT-Ki-Q^+irf,.      (7.35) 

/  /  (  / 

Using  (7.29),  we  rewrite  this  expression  as 

^\d[C,jdj   -   d[C,jd].  (7.36) 

In  (7.36),  we  reverse  the  roles  of  i  and  j  and  take  the  average  of  the  resulting  expression 
and  (7.36),  obtaining 


-j-Zid^-d^yc.jid^-dj)  s  0. 


(7.37) 


This  completes  the  proof  of  proposition  7.1. 


As  a  very  simple  special  case,  we  consider  an  immersed  circle  under  tension.  That 
is,  the  initial  configuration  of  the  boundary  is  given  by 


U')   = 


Acos— 
a 

i4sm— 


(7.38) 
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where  A>fl>0  and  j€[0;2iTa].  We  show  that  a  solution  of  (7.7)  is  given  by 


u  =  Q 


P  = 


5(1- X)  if  llill^A 
a     A 
0      otherwise 


i(j,r)  =  2Co(') 
At  time  0,  the  tension  T  is 


Since 


r  =  S(— -1). 
a 


'-sin(f)' 

cos(f) 
a 


(7.39) 
(7.40) 

(7.41) 
(7.42) 


(7.43) 


the  distribution  £  at  time  0  is  given  by 

2ira 

'( 

Q       a         a  a  a  a  a 

We  show  that  this  is  the  distributional  gradient  of  the  function  p  defined  by  (7.40). 
By  definition, 


ATra 

(]E.,m)=  /5(— -l)-(-cos(-),~sin(-)ii:(Acos-^sin-)d:f.  (7.44) 


(S>,m:)  :=  -ipS-yn.)- 
Using  (7.40)  and  the  divergence  theorem,  we  obtain 

(SP.1K)  =  -   /  5(i-i)£ii:(i)4i 


(7.45) 


i 


M^ 


a     A 


(Sp,ii:)=  -5(^-4")    f  aarft. 

Ikll-A 


a      A 


(7.46) 

(7.47) 

(7.47)  is  equal  to  (7.44). 

In  summary,  an  immersed  drde  under  tension  rests  in  the  fluid,  it  is  prevented 
from  shrinking  by  an  increased  pressure  in  the  interior  of  the  circle  and  the  pressure 

difference  is  5(— -  -7-). 
a      A 
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The  drclc  does  not  shrink  because  of  the  incompressibility.  Therefore,  this  problem 
provides  an  opportunity  to  test  how  well  the  incompressibility  condition  is  approximated 
by  the  numerical  method.  If  the  drcle  shrinks  in  the  numerical  computations,  there  are 
two  possible  explanations.  Either  markers  which  are  inside  the  drcle  at  time  0  leave  the 
circle  -  there  is  nothing  in  the  method  that  explicitly  prevents  this  -,  or  the  markers  in 
the  circle  are  more  densely  packed  at  later  times  than  at  time  0,  in  violation  of  the  in- 
compressibility condition.  We  have  found  that  shrinking  never  occurs  for  the  first  rea- 
son. Markers  which  are  inside  the  drcle  at  time  zero  remain  inside  the  drde.  Shrinking 
does  occur  for  the  second  reason.  The  following  figure  shows  the  results  of  a  computa- 
tion in  which  the  initial  configuiation  was  a  drde  around  (0.5,0.5)  with  a  radius 
0.30625.  iks  is  chosen  such  that  the  drde  with  radius  0.25  is  a  resting  configuration  of 
the  boundary.  The  curves  shew  the  area  induded  by  the  computed  boundary  plotted 
against  time  for  different  numbers  N  of  markers.  The  lowest  curve  corresponds  to 
N=100,M=36,  the  middle  one  to  N=400,M^76,  the  upper  one  to  N=  1600,M=  156. 

Ar=-f-,  v=1.0  and  S=0.5.    We  use  the  modified  fractional  step  method  (6.1b)  with 

two  iterations  of  the  projection  algorithm  in  each  step.  All  scalar  Helmholtz  and  Poisson 
fjToblems  are  solved  up  to  rounding  errors;  sec  the  remark  below  concerning  this  point. 
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The  amoimts  of  area  lost  at  time  1  are  0.00150  for  N=100,  0.00094  for  N=400  and 
0.00045  for  N=  1600.  Thus  the  loss  of  area  appears  to  decrease  for  increasing  N. 

It  is  crucial  to  use  more  than  one  iteration  for  each  Helmholtz  or  Poisson  problem, 
since  the  mesh  is  quite  degenerate  along  the  immersed  boundary.  We  illustrate  this  by 
the  following  figure,  showing  the  Voronoi  diagram  at  time  0  for  N=100.  This  degen- 
eracy causes  some  problems  for  our  two-level  algorithm  because  it  uses  a  uniform  auzili- 

.     f  <  * 

ary  grid.  One  might,  of  course,  consider  using  auxiliary  grids  which  also  have  a  higher 
concentration  of  points  along  the  immersed  boundary  than  elsewhere.  We  have  not  yet 
tested  this  possibility. 


For  large  S,  we  have  observed  unstable  behavior.  After  shrinking  slowly  for  a  few 
time  steps,  the  circle  suddenly  breaks,  rapidly  shrinking  down  to  the  resting  radius  of 
0.25. 
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8.  Open  problems  and  plans  for  ftatnre  work 

The  weakest  piece  m  the  method  which  we  have  studied  is  probably  the  projection 
algoritfain  5.1.  The  fact  that  Ihe  algorithm,  when  run  over  several  iterations,  sometimes 
requires  a  large  power  k  in  the  smoothing  factor 

(/  -t-  tiiiy, 
indicates  a  stability  problem  wliidi  is  not  satisfactorily  removed  by  the  smoothing  intro- 
duced in  section  5.  At  the  present  stage  of  this  work,  the  smoothing  step  is  useful 
mainly  as  a  guarantee  that  n;"i;ung  v.^J!  gc  wrong  within  the  first  and  second  iteration  of 
the  projection  algorithm.  It  do  *  net  really  make  it  feasible  to  run  the  algorithm  to  con- 
vergence. To  overcome  thi;}  'Ji^rJculty,  the  definitions  of  the  operators  D  and  Q  should 
perhaps  be  modified.  The  experiment  with  Dj  presented  in  section  5,  and  some  similar 
experiments  about  which  wc  h^^vc  not  reported  here,  indicate,  however,  that  slight 
changes  in  the  definitions  of  D  and  i2  aiay  lead  to  a  sharp  decrease  in  accuracy  in  the 
projection  step,  possibly  even  to  a  complete  loss  of  pointwise  convergence.  For  a  better 
understanding  of  this  observation,  a  proof  of  the  convergence  shown  in  table  5.1  may  be 
crucial. 

In  aU  the  experiments  mentioned  in  sections  6  and  7,  we  have  taken  At  to  be  5/iv, 
where  8  was  a  constant  not  much  smaller  than  1.  We  have  observed  cases  in  which 
reduction  of  8  increases  the  error,  in  a  way  similar  to  the  increase  in  error  observed 
when  decreasing  the  viscosity  coefficient  v.  This  may  have  the  following  explanation. 
The  projection  step  may  amplify  certain  high  frequency  errors  which  do  not  always  des- 
troy the  accuracy  because  they  are  sufficiently  damped  by  the  smoothing  provided  by  the 
viscosity.  In  the  fractional  step  method  (6.1),  this  smoothing  appears  in  the  form 

i"*i  :=  (/-uZ,")- V  +  Arf  *i), 
where  (i)  =  vAr.  If  N  is  fixed  and  u)  is  decreased,  the  smoothing  step  becomes  less  effec- 
tive and  may  fail  to  damp  the  high  frequency  errors  introduced  by  the  projection  step 
sufficiently. 
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The  accuracy  of  the  method  should  be  increased.  One  can  easily  and  in  various 
different  ways  obtain  0(h^)  in  (3.41)  for  the  divergence  and  gradient  operators,  thus 
obtaining  pointwise  consistent  discretizations.  It  is,  of  course,  not  certain  that  such 
operators  are  superior  to  the  weakly  first  order  consistent  operators  D  ,G.  which  we  have 
used.  Vf'e  have  tested  one  such  possibility.  It  turned  out  to  give  worse  results  than  D,Q 
when  used  within  the  projection  algorithm  in  section  5.  Pointwise  consistent  discretiza- 
tions of  more  than  first  order  woidd  probably  require  larger  stencils. 

The  multigrid  method  described  in  section  4  has  the  disadvantage  of  not  being  vec- 
torizable.  It  may  therefore  be  useful  to  replace  i\t  relaxation  scheme  by  an  analogous 
scheme  using  several  colors.  On  the  torus,  the  minimum  possible  number  of  colors  is  7; 
see  Hilbert  and  Cohn-Vossen  (1932).  The  proof  that  7  colors  suffice  easily  leads  to  an 
efficient  algorithm  for  the  computation  of  such  colorings.  It  may,  however,  be  sufficient 
to  use  a  partitioning  into  less  than  7  colors  so  that  .the  coupling  in  the  discrete  Helmhcltz 
operator  between  points  belonging  to  one  color  is  not  niscessarily  zero  but  relatively 
weak.  A  relaxation  sweep  through  one  such  color  could  then  be  carried  out  Jacobi-like, 
and  could  be  completely  parallelized.  A  similar  algorithm  has  been  used  succes'^fully  to 
solve  certain  difference  equations  on  regular  grids,  using  a  compact  9-point  stencil;  see 
Stiiben  and  Trottenberg  (1981),  p.  113.  In  this  algorithm,  only  two  different  colors  are 
used,  with  weak  coupling  within  each  color. 

We  also  plan  further  work  on  the  algorithm  discussed  in  section  7,  an  implementa- 
tion of  the  method  with  the  boundary  condition 

u  =  0  on  dft 
and  an  implementation  in  three  space  dimensions. 
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